#334 compute geometry bbox in target crs if possible

Author: jdries

openeo_driver/dry_run.py

@@ -39,6 +39,7 @@
 
 import numpy
 import shapely.geometry.base
+from pyproj import CRS
 from shapely.geometry import Point, Polygon, MultiPolygon, GeometryCollection
 from shapely.geometry.base import BaseGeometry
 
@@ -53,7 +54,7 @@
     AggregatePolygonResult,
     AggregatePolygonSpatialResult,
 )
-from openeo_driver.util.geometry import geojson_to_geometry, GeometryBufferer
+from openeo_driver.util.geometry import geojson_to_geometry, GeometryBufferer, BoundingBox
 from openeo_driver.utils import to_hashable, EvalEnv
 
 _log = logging.getLogger(__name__)
@@ -519,7 +520,11 @@ def filter_bbox(
                                          "crs": (crs or "EPSG:4326")})
 
     def filter_spatial(self, geometries):
-        geometries, bbox = self._normalize_geometry(geometries)
+        crs = 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)
         cube = self.filter_bbox(**bbox, operation="weak_spatial_extent")
         return cube._process(operation="filter_spatial", arguments={"geometries": geometries})
 
@@ -570,11 +575,15 @@ def aggregate_spatial(
         geoms_is_empty = isinstance(geometries, DriverVectorCube) and len(geometries.get_geometries()) == 0
         cube = self
         if not geoms_is_empty:
-            geometries, bbox = self._normalize_geometry(geometries)
+            crs = 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)
             cube = self.filter_bbox(**bbox, operation="weak_spatial_extent")
         return cube._process(operation="aggregate_spatial", arguments={"geometries": geometries})
 
-    def _normalize_geometry(self, geometries) -> Tuple[Union[DriverVectorCube, DelayedVector, BaseGeometry], dict]:
+    def _normalize_geometry(self, geometries, target_crs=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)
@@ -586,8 +595,17 @@ def _normalize_geometry(self, geometries) -> Tuple[Union[DriverVectorCube, Delay
             # TODO: buffer distance of 10m assumes certain resolution (e.g. sentinel2 pixels)
             # TODO: use proper distance for collection resolution instead of using a default distance?
             # TODO: or eliminate need for buffering in the first place? https://github.com/Open-EO/openeo-python-driver/issues/148
-            bbox = geometries.buffer_points(distance=10).get_bounding_box()
-            crs = geometries.get_crs_str()
+            if target_crs is not None:
+                is_utm = crs == "AUTO:42001" or "Auto42001" in str(target_crs)
+                if is_utm:
+                    target_crs = BoundingBox.from_wsen_tuple(geometries.get_bounding_box(),crs=geometries.get_crs()).best_utm()
+                else:
+                    target_crs = BoundingBox.normalize_crs(target_crs)
+                bbox = geometries.buffer_points(distance=10).reproject(CRS.from_user_input( target_crs )).get_bounding_box()
+                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))
         elif isinstance(geometries, str):

tests/test_dry_run.py

@@ -2185,3 +2185,57 @@ def test_complex_extract_load_stac(dry_run_env,dry_run_tracer):
     assert (loadparams.global_extent == expected_extent)
     assert (loadparams.target_resolution == [10,10])
 
+def test_normalize_geometries(dry_run_env,dry_run_tracer):
+    geometry_path = str(get_path("geojson/points.geojson"))
+    pg = {
+        "lc": {"process_id": "load_collection", "arguments": {"id": "S2_FOOBAR"}},
+        "vector": {"process_id": "read_vector", "arguments": {"filename": geometry_path}},
+        "agg": {
+            "process_id": "aggregate_spatial",
+            "arguments": {
+                "data": {"from_node": "lc"},
+                "geometries": {"from_node": "vector"},
+                "reducer": {
+                    "process_graph": {
+                        "mean": {
+                            "process_id": "mean",
+                            "arguments": {"data": {"from_parameter": "data"}},
+                            "result": True,
+                        }
+                    }
+                },
+            },
+        },
+        "raster": {
+            "process_id": "vector_to_raster",
+            "arguments": {"data": {"from_node": "agg"}, "target_data_cube": {"from_node": "lc"}},
+        },
+        "rename": {
+            "process_id": "rename_labels",
+            "arguments": {"data": {"from_node": "raster"}, "dimension": "bands", "target": ["score"]},
+            "result": True,
+        },
+    }
+    save_result = evaluate(pg, env=dry_run_env)
+    source_constraints = dry_run_tracer.get_source_constraints(merge=True)
+    print(source_constraints)
+
+    dry_run_env = dry_run_env.push({ENV_SOURCE_CONSTRAINTS: source_constraints})
+    params = [_extract_load_parameters(dry_run_env,source_id) for source_id, _ in source_constraints]
+    extents = [BoundingBox.from_dict(param.spatial_extent).as_polygon() for param in params]
+
+    from shapely.geometry import mapping
+    from shapely.ops import transform
+    import json
+    import pyproj
+    t = pyproj.Transformer.from_crs(pyproj.CRS.from_user_input(32632), pyproj.CRS.from_user_input(4326), always_xy=True)
+    extents += [transform(t.transform,BoundingBox.from_dict(param.global_extent).as_polygon().segmentize(1000)) for param in params]
+
+    #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}