aggregate_spatial with points: more robust buffering approach

Author: jdries

openeo_driver/dry_run.py

@@ -561,10 +561,12 @@ def filter_bbox(
 
     def filter_spatial(self, geometries):
         crs = None
+        resolution = None
         if len(self.metadata.spatial_dimensions) > 0:
             spatial_dim = self.metadata.spatial_dimensions[0]
             crs = spatial_dim.crs
-        geometries, bbox = self._normalize_geometry(geometries, target_crs=crs)
+            resolution = spatial_dim.step
+        geometries, bbox = self._normalize_geometry(geometries, target_crs=crs, target_resolution=resolution)
         cube = self.filter_bbox(**bbox, operation="weak_spatial_extent")
         return cube._process(operation="filter_spatial", arguments={"geometries": geometries})
 
@@ -622,19 +624,25 @@ def aggregate_spatial(
         cube = self
         if not geoms_is_empty:
             crs = None
+            resolution = None
             if len(self.metadata.spatial_dimensions) > 0:
                 spatial_dim = self.metadata.spatial_dimensions[0]
+                resolution = spatial_dim.step
                 crs = spatial_dim.crs
-            geometries, bbox = self._normalize_geometry(geometries, target_crs=crs)
+            geometries, bbox = self._normalize_geometry(geometries, target_crs=crs, target_resolution=resolution)
             cube = self.filter_bbox(**bbox, operation="weak_spatial_extent")
         return cube._process(operation="aggregate_spatial", arguments={"geometries": geometries})
 
     def _normalize_geometry(
-        self, geometries, target_crs=None
+        self, geometries, target_crs=None, target_resolution = None
     ) -> Tuple[Union[DriverVectorCube, DelayedVector, BaseGeometry], dict]:
         """
         Helper to preprocess geometries (as used in aggregate_spatial and mask_polygon)
         and extract bbox (e.g. for filter_bbox)
+
+        :param geometries: geometries as BaseGeometry, GeoJSON dict, DelayedVector or DriverVectorCube
+        :param target_crs: target CRS to reproject geometries to
+        :param target_resolution: target resolution for geometries, in units of the target CRS, None by default which will use 10m
         """
         _log.debug(f"_normalize_geometry with {type(geometries)}")
         # TODO #71 #114 EP-3981 normalize to vector cube instead of GeometryCollection
@@ -660,21 +668,25 @@ def _normalize_geometry(
 
                 if len(points) > 0:
                     point_hull = reproject_geometry(shapely.geometry.MultiPoint(points).envelope,geometries.get_crs(),CRS.from_user_input(target_crs))
-                    loi_point = point_hull.representative_point()
-                    bufferer = GeometryBufferer.from_meter_for_crs(
-                        distance=10, crs=target_crs, loi=(loi_point.x,loi_point.y), loi_crs=target_crs
-                    )
-                    buffered_points = bufferer.buffer(point_hull)
+                    buffered_points = None
+                    if target_resolution == None:
+                        loi_point = point_hull.representative_point()
+                        bufferer = GeometryBufferer.from_meter_for_crs(
+                            distance=10, crs=target_crs, loi=(loi_point.x,loi_point.y), loi_crs=target_crs
+                        )
+                        buffered_points = bufferer.buffer(point_hull)
+                    else:
+                        buffered_points = point_hull.buffer(distance=target_resolution)
                     other_hull = other_hull.union(buffered_points)
                 bbox = other_hull.bounds
                 crs = target_crs
             else:
                 bbox = geometries.buffer_points(distance=10).get_bounding_box()
                 crs = geometries.get_crs_str()
         elif isinstance(geometries, dict):
-            return self._normalize_geometry(geojson_to_geometry(geometries))
+            return self._normalize_geometry(geojson_to_geometry(geometries),target_crs, target_resolution)
         elif isinstance(geometries, str):
-            return self._normalize_geometry(DelayedVector(geometries))
+            return self._normalize_geometry(DelayedVector(geometries),target_crs, target_resolution)
         elif isinstance(geometries, DelayedVector):
             bbox = geometries.bounds
         elif isinstance(geometries, shapely.geometry.base.BaseGeometry):

tests/test_dry_run.py

@@ -652,10 +652,10 @@ def test_aggregate_spatial_only(dry_run_env, dry_run_tracer):
             approxify(
                 {
                     "crs": "EPSG:32631",
-                    "east": 610869.8770073673,
-                    "north": 552758.5425218772,
-                    "south": 331634.0169357686,
-                    "west": 388860.8935124034,
+                    "east": 611139.132350062,
+                    "north": 552758.6209151235,
+                    "south": 331633.9380733739,
+                    "west": 388860.86764994,
                 }
             ),
         ),
@@ -2376,7 +2376,7 @@ def test_complex_extract_load_stac(dry_run_env, dry_run_tracer):
     loadparams = _extract_load_parameters(dry_run_env, source_id_bands)
 
     print(loadparams)
-    expected_extent = {"crs": "EPSG:32633", "east": 400380, "north": 4700400, "south": 4679610, "west": 380120}
+    expected_extent = {'west': 379840, 'south': 4679580, 'east': 400460, 'north': 4700510, 'crs': 'EPSG:32633'}
     assert loadparams.global_extent == expected_extent
     assert loadparams.bands == ["B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B11", "B12"]
     assert loadparams.pixel_buffer == None
@@ -2441,12 +2441,11 @@ def test_normalize_geometries(dry_run_env, dry_run_tracer):
 
     # print(json.dumps(dict(type="FeatureCollection",features=[dict(type="Feature",geometry=mapping(e),properties={}) for e in extents])))
 
-    assert params[0].global_extent == {
-        "crs": "EPSG:32632",
-        "east": 636980,
-        "north": 5729350,
-        "south": 5654910,
-        "west": 89920,
+    assert params[0].global_extent == {'crs': 'EPSG:32632',
+     'east': 638890,
+     'north': 5736520,
+     'south': 5642750,
+     'west': 84850
     }