FIX: Force 2D geometries in rasters.crop as odc.geo.xr.crop don't…

… work with 3D geometries

CHANGES.md

@@ -2,6 +2,7 @@
 
 ## 1.xx.x 
 
+- FIX: Force 2D geometries in `rasters.crop` as `odc.geo.xr.crop` don't work with 3D geometries
 - OPTIM: Only write rasters on disk with `windowed=True` in case of big rasters _(w/*h > 20,000*20,000 pixels)_
 - CI: Refactor `.gitlab-ci.yml` file with new GitLab templates
 

ci/test_geometry.py

@@ -192,3 +192,19 @@ def test_nearest_neighbors():
     closest, distances = nearest_neighbors(
         src, candidates, method="k_neighbors", k_neighbors=len(candidates) + 10
     )
+
+
+def test_force_2_or_3d():
+    """Force 2D or 3D"""
+    aoi_kml = vectors.read(vectors_path().joinpath("aoi.kml"))
+
+    assert aoi_kml.has_z.any()
+
+    # Check equality
+    ci.assert_geom_equal(
+        aoi_kml, geometry.force_3d(geometry.force_2d(aoi_kml)), ignore_z=False
+    )
+
+    # Check Z type
+    assert not geometry.force_2d(aoi_kml).has_z.any()
+    assert geometry.force_3d(geometry.force_2d(aoi_kml)).has_z.any()

ci/test_rasters.py

@@ -26,7 +26,7 @@
 import xarray as xr
 
 from ci.script_utils import KAPUT_KWARGS, dask_env, get_output, rasters_path, s3_env
-from sertit import ci, path, rasters, unistra, vectors
+from sertit import ci, geometry, path, rasters, unistra, vectors
 from sertit.rasters import (
     FLOAT_NODATA,
     INT8_NODATA,
@@ -330,6 +330,13 @@ def test_crop(tmp_path, xda, xds, xda_dask, mask):
     ci.assert_raster_equal(xda_cropped, raster_cropped_xarray_path)
     ci.assert_raster_equal(xds_cropped, raster_cropped_xarray_path)
 
+    # Test with mask with Z
+    mask_z = geometry.force_3d(mask)
+    crop_z = rasters.crop(xda, mask_z)
+    assert crop_z.chunks is not None
+    np.testing.assert_array_equal(crop_xda, crop_z)
+    ci.assert_xr_encoding_attrs(crop_xda, crop_z)
+
 
 @s3_env
 @dask_env

sertit/geometry.py

@@ -603,3 +603,21 @@ def _get_radius_nearest(
 
     # Return indices and distances
     return closest, closest_dist
+
+
+def force_2d(vect: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
+    """Have the force_2d function even with geopandas < 1.0.0"""
+    try:
+        return vect.force_2d()
+    except AttributeError:
+        vect.geometry = vect.geometry.apply(shapely.force_2d)
+        return vect
+
+
+def force_3d(vect: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
+    """Have the force_3d function even with geopandas < 1.0.0"""
+    try:
+        return vect.force_3d()
+    except AttributeError:
+        vect.geometry = vect.geometry.apply(shapely.force_3d)
+        return vect

sertit/rasters.py

@@ -771,7 +771,7 @@ def crop(
 
         # Convert the shapes in the right format
         if isinstance(shapes, (gpd.GeoDataFrame, gpd.GeoSeries)):
-            shapes = shapes.to_crs(xds.rio.crs)
+            shapes = geometry.force_2d(shapes.to_crs(xds.rio.crs))
 
             try:
                 shapes = shapes.union_all()