Repository: corteva/rioxarray

Issue:
In to_raster, check for original raster object profile for writing to raster
![image](https://user-images.githubusercontent.com/8699967/60894486-62b60b80-a228-11e9-93a2-eb11ff2580a7.png)


Discussion:

Patch:
diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py
index acb5a1f4..1d791fc8 100644
--- a/rioxarray/rioxarray.py
+++ b/rioxarray/rioxarray.py
@@ -375,11 +375,16 @@ class RasterArray(XRasterBase):
                     self._obj.attrs.get("fill_value", self._obj.attrs.get("nodata")),
                 ),
             )
+
+        # look in places used by `xarray.open_rasterio`
         if self._nodata is None:
             try:
-                self._nodata = self._obj.attrs["nodatavals"][0]
-            except (KeyError, IndexError):
-                pass
+                self._nodata = self._obj._file_obj.acquire().nodata
+            except AttributeError:
+                try:
+                    self._nodata = self._obj.attrs["nodatavals"][0]
+                except (KeyError, IndexError):
+                    pass
 
         if self._nodata is None:
             self._nodata = False
@@ -750,7 +755,7 @@ class RasterArray(XRasterBase):
 
         return cl_array
 
-    def clip(self, geometries, crs, all_touched=False):
+    def clip(self, geometries, crs, all_touched=False, drop=True):
         """
         Crops a :class:`xarray.DataArray` by geojson like geometry dicts.
 
@@ -762,10 +767,14 @@ class RasterArray(XRasterBase):
             A list of geojson geometry dicts.
         crs: :obj:`rasterio.crs.CRS`
             The CRS of the input geometries.
-        all_touched : boolean, optional
+        all_touched : bool, optional
             If True, all pixels touched by geometries will be burned in.  If
             false, only pixels whose center is within the polygon or that
             are selected by Bresenham's line algorithm will be burned in.
+        drop: bool, optional
+            If True, drop the data outside of the extent of the mask geoemtries
+            Otherwise, it will return the same raster with the data masked.
+            Default is True.
 
         Returns
         -------
@@ -786,7 +795,7 @@ class RasterArray(XRasterBase):
             >>> cropped = xds.rio.clip(geometries=cropping_geometries, crs=4326)
         """
         geometries = [
-            rasterio.warp.transform_geom(self.crs, CRS.from_user_input(crs), geometry)
+            rasterio.warp.transform_geom(CRS.from_user_input(crs), self.crs, geometry)
             for geometry in geometries
         ]
 
@@ -806,8 +815,11 @@ class RasterArray(XRasterBase):
             },
             dims=(self.y_dim, self.x_dim),
         )
+        cropped_ds = self._obj.where(clip_mask_xray, drop=drop)
+        if self.nodata is not None and not np.isnan(self.nodata):
+            cropped_ds = cropped_ds.fillna(self.nodata)
 
-        cropped_ds = self._obj.where(clip_mask_xray, drop=True).astype(self._obj.dtype)
+        cropped_ds = cropped_ds.astype(self._obj.dtype)
 
         if (
             cropped_ds.coords[self.x_dim].size < 1
@@ -918,7 +930,7 @@ class RasterArray(XRasterBase):
         **profile_kwargs
             Additional keyword arguments to pass into writing the raster. The
             nodata, transform, crs, count, width, and height attributes
-            are automatically added.
+            are ignored.
 
         """
         width, height = self.shape
@@ -927,6 +939,32 @@ class RasterArray(XRasterBase):
         count = 1
         if extra_dim is not None:
             count = self._obj[extra_dim].size
+
+        # get the output profile from the rasterio object
+        # if opened with xarray.open_rasterio()
+        try:
+            out_profile = self._obj._file_obj.acquire().profile
+        except AttributeError:
+            out_profile = {}
+        out_profile.update(profile_kwargs)
+
+        # filter out the generated attributes
+        out_profile = {
+            key: value
+            for key, value in out_profile.items()
+            if key
+            not in (
+                "driver",
+                "height",
+                "width",
+                "crs",
+                "transform",
+                "nodata",
+                "count",
+                "dtype",
+            )
+        }
+
         with rasterio.open(
             raster_path,
             "w",
@@ -940,7 +978,7 @@ class RasterArray(XRasterBase):
             nodata=(
                 self.encoded_nodata if self.encoded_nodata is not None else self.nodata
             ),
-            **profile_kwargs,
+            **out_profile,
         ) as dst:
 
             if self.encoded_nodata is not None:
@@ -1092,7 +1130,7 @@ class RasterDataset(XRasterBase):
         clipped_dataset.attrs["creation_date"] = str(datetime.utcnow())
         return clipped_dataset
 
-    def clip(self, geometries, crs, all_touched=False):
+    def clip(self, geometries, crs, all_touched=False, drop=True):
         """
         Crops a :class:`xarray.Dataset` by geojson like geometry dicts.
 
@@ -1111,6 +1149,10 @@ class RasterDataset(XRasterBase):
             If True, all pixels touched by geometries will be burned in.  If
             false, only pixels whose center is within the polygon or that
             are selected by Bresenham's line algorithm will be burned in.
+        drop: bool, optional
+            If True, drop the data outside of the extent of the mask geoemtries
+            Otherwise, it will return the same raster with the data masked.
+            Default is True.
 
         Returns
         -------
@@ -1133,7 +1175,7 @@ class RasterDataset(XRasterBase):
         clipped_dataset = xarray.Dataset(attrs=self._obj.attrs)
         for var in self.vars:
             clipped_dataset[var] = self._obj[var].rio.clip(
-                geometries, crs=crs, all_touched=all_touched
+                geometries, crs=crs, all_touched=all_touched, drop=drop
             )
         clipped_dataset.attrs["creation_date"] = str(datetime.utcnow())
         return clipped_dataset
diff --git a/sphinx/history.rst b/sphinx/history.rst
index 1ce552aa..6c761b32 100644
--- a/sphinx/history.rst
+++ b/sphinx/history.rst
@@ -5,6 +5,10 @@ History
 -----
 - Add support for scalar coordinates in reproject (issue #15)
 - Updated writing encoding for FutureWarning (issue #18)
+- Use input raster profile for defaults to write output raster profile if opened with `xarray.open_rasterio` (issue #19)
+- Preserve None nodata if opened with `xarray.open_rasterio` (issue #20)
+- Added `drop` argument for `clip()` (issue #25)
+- Fix order of `CRS` for reprojecting geometries in `clip()` (pull #24)
 
 0.0.5
 -----


Test patch:
diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py
index 8c479c3a..049b15a3 100644
--- a/test/integration/test_integration_rioxarray.py
+++ b/test/integration/test_integration_rioxarray.py
@@ -4,6 +4,7 @@ import numpy
 import pytest
 import rasterio
 import xarray
+from affine import Affine
 from numpy.testing import assert_almost_equal, assert_array_equal
 from rasterio.crs import CRS
 
@@ -261,10 +262,9 @@ def test_clip_geojson():
         # add transform for test
         comp_subset.attrs["transform"] = tuple(comp_subset.rio.transform(recalc=True))
         # add grid mapping for test
-        comp_subset.attrs["grid_mapping"] = "spatial_ref"
-        comp_subset.coords["spatial_ref"] = 0
+        comp_subset.rio.write_crs(subset.rio.crs, inplace=True)
         # make sure nodata exists for test
-        comp_subset.attrs["_FillValue"] = comp_subset.attrs["nodatavals"][0]
+        comp_subset.attrs["_FillValue"] = comp_subset.rio.nodata
 
         geometries = [
             {
@@ -282,7 +282,7 @@ def test_clip_geojson():
         ]
 
         # test data array
-        clipped = xdi.rio.clip(geometries, subset.rio.crs)
+        clipped = xdi.rio.clip(geometries, comp_subset.rio.crs)
         _assert_xarrays_equal(clipped, comp_subset)
 
         # test dataset
@@ -293,6 +293,39 @@ def test_clip_geojson():
         _assert_xarrays_equal(clipped_ds, comp_subset_ds)
 
 
+def test_clip_geojson__no_drop():
+    with xarray.open_rasterio(
+        os.path.join(TEST_COMPARE_DATA_DIR, "small_dem_3m_merged.tif")
+    ) as xdi:
+        geometries = [
+            {
+                "type": "Polygon",
+                "coordinates": [
+                    [
+                        [-93.880889448126, 41.68465068553298],
+                        [-93.89966980835203, 41.68465068553298],
+                        [-93.89966980835203, 41.689430423525266],
+                        [-93.880889448126, 41.689430423525266],
+                        [-93.880889448126, 41.68465068553298],
+                    ]
+                ],
+            }
+        ]
+        # test data array
+        clipped = xdi.rio.clip(geometries, "epsg:4326", drop=False)
+        assert clipped.rio.crs == xdi.rio.crs
+        assert clipped.shape == xdi.shape
+        assert clipped.sum().item() == 2150801411
+
+        # test dataset
+        clipped_ds = xdi.to_dataset(name="test_data").rio.clip(
+            geometries, "epsg:4326", drop=False
+        )
+        assert clipped_ds.rio.crs == xdi.rio.crs
+        assert clipped_ds.test_data.shape == xdi.shape
+        assert clipped_ds.test_data.sum().item() == 2150801411
+
+
 def test_transform_bounds():
     with xarray.open_dataarray(
         os.path.join(TEST_INPUT_DATA_DIR, "MODIS_ARRAY.nc")
@@ -729,6 +762,40 @@ def test_to_raster_3d(tmpdir):
         assert_array_equal(rds.read(), xds.values)
 
 
+def test_to_raster__preserve_profile__none_nodata(tmpdir):
+    tmp_raster = tmpdir.join("output_profile.tif")
+    input_raster = tmpdir.join("input_profile.tif")
+
+    transform = Affine.from_gdal(0, 512, 0, 0, 0, 512)
+    with rasterio.open(
+        str(input_raster),
+        "w",
+        driver="GTiff",
+        height=512,
+        width=512,
+        count=1,
+        crs="+init=epsg:4326",
+        transform=transform,
+        dtype=rasterio.float32,
+        tiled=True,
+        tilexsize=256,
+        tileysize=256,
+    ) as rds:
+        rds.write(numpy.empty((1, 512, 512), dtype=numpy.float32))
+
+    with xarray.open_rasterio(str(input_raster)) as mda:
+        mda.rio.to_raster(str(tmp_raster))
+
+    with rasterio.open(str(tmp_raster)) as rds, rasterio.open(str(input_raster)) as rdc:
+        assert rds.count == rdc.count
+        assert rds.crs == rdc.crs
+        assert_array_equal(rds.transform, rdc.transform)
+        assert_array_equal(rds.nodata, rdc.nodata)
+        assert_array_equal(rds.read(), rdc.read())
+        assert rds.profile == rdc.profile
+        assert rds.nodata is None
+
+
 def test_missing_dimensions():
     test_ds = xarray.Dataset()
     with pytest.raises(KeyError):