Repository: corteva/rioxarray

Issue:
add `drop` option to `clip()`
Looks like the clip function in QGIS defaults to match the extent of the masked geometry:
![image](https://user-images.githubusercontent.com/8699967/61016683-1a950700-a356-11e9-8606-fff8b82d8e44.png)

But, it also has the option to keep the original resolution of the original raster. I am thinking about adding the `drop` argument to the `clip()` function to preserve the original raster an only return a masked version.

Related to ihttps://github.com/pydata/xarray/issues/501

Discussion:
snowman2: Thanks @alfredoahds!

Patch:
diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py
index a067afc4..1d791fc8 100644
--- a/rioxarray/rioxarray.py
+++ b/rioxarray/rioxarray.py
@@ -755,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.
 
@@ -767,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
         -------
@@ -791,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
         ]
 
@@ -811,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
@@ -1123,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.
 
@@ -1142,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
         -------
@@ -1164,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 5dfe9aed..6c761b32 100644
--- a/sphinx/history.rst
+++ b/sphinx/history.rst
@@ -7,6 +7,8 @@ History
 - 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 857c23d9..049b15a3 100644
--- a/test/integration/test_integration_rioxarray.py
+++ b/test/integration/test_integration_rioxarray.py
@@ -262,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 = [
             {
@@ -283,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
@@ -294,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")