added 'drop' argument for clip; fixed geometry crs reprojection in clip

rioxarray/rioxarray.py

@@ -755,7 +755,7 @@ def clip_box(self, minx, miny, maxx, maxy, auto_expand=False, auto_expand_limit=
 
         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 @@ def clip(self, geometries, crs, all_touched=False):
             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 @@ def clip(self, geometries, crs, all_touched=False):
             >>> 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 @@ def clip(self, geometries, crs, all_touched=False):
             },
             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 @@ def clip_box(self, minx, miny, maxx, maxy, auto_expand=False, auto_expand_limit=
         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 @@ def clip(self, geometries, crs, all_touched=False):
             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 @@ def clip(self, geometries, crs, all_touched=False):
         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

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/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")