1 files changed, 41 insertions, 19 deletions

diff --git a/common_gdal.py b/common_gdal.py
index 4fea54a..35511a3 100644
--- a/common_gdal.py
+++ b/common_gdal.py
@@ -94,15 +94,11 @@ def getSRS(srs_str : Optional[str]) -> osr.SpatialReference:
     return srs
 
 def getExtent(extent : Optional[tuple[float, float, float, float]],
-              srs : Optional[osr.SpatialReference] = None) -> tuple[ogr.Geometry, ogr.Geometry]:
+              srs : Optional[osr.SpatialReference] = None) -> ogr.Geometry:
     """Convert extent (minX, minY, maxX, maxY) into a polygon and assign the
-    given SRS.  Return a pair with the densified and non-densified extent.
-    Like apps/ogr2ogr_lib.cpp, the former is obtained by segmentizing the
-    polygon to make sure it is sufficiently densified when transforming to
-    source layer SRS for spatial filtering."""
-
+    given SRS."""
     if extent is None:
-        return None, None
+        return None
 
     if not isinstance(extent, tuple) or len(extent) != 4:
         raise RuntimeError(f'Invalid extent {extent}')
@@ -130,18 +126,44 @@ def getExtent(extent : Optional[tuple[float, float, float, float]],
     if not srs2.IsSame(srs):
         polygon.TransformTo(srs)
 
-    # densify the rectangle to avoid issues when reprojecting to the
-    # source layer SRS, cf. apps/ogr2ogr_lib.cpp:ApplySpatialFilter()
-    polygon_dense = polygon.Clone()
-    segment_distance_metre = 10 * 1000
-    if srs.IsGeographic():
-        dfMaxLength = segment_distance_metre / math.radians(srs.GetSemiMajor())
-        polygon_dense.Segmentize(dfMaxLength)
-    elif srs.IsProjected():
-        dfMaxLength = segment_distance_metre / srs.GetLinearUnits()
-        polygon_dense.Segmentize(dfMaxLength)
-
-    return polygon_dense, polygon
+    return polygon
+
+def getSpatialFilterFromGeometry(geom : ogr.Geometry,
+                                 srs : osr.SpatialReference,
+                                 fail_if_srs_not_equivalent : bool = False) -> ogr.Geometry:
+    """Make the geometry suitable for use as a spatial filter.  It is
+    densified the SRS:s are not equivalent, unless fail_if_srs_not_equivalent is True in
+    which case it errors out (this ensure clipping geometries remain simple)."""
+    cloned = False
+    geom_srs = geom.GetSpatialReference()
+    if not geom_srs.IsSame(srs, [ 'IGNORE_DATA_AXIS_TO_SRS_AXIS_MAPPING=YES',
+                                  'CRITERION=EQUIVALENT']):
+        if fail_if_srs_not_equivalent:
+            raise RuntimeError(f'Geometry {geom.ExportToWkt()} has SRS ' +
+                                 geom_srs.ExportToPrettyWkt() + ',\nexpected ' +
+                                 srs.ExportToPrettyWkt())
+        # densify the geometry (a rectangle) to avoid issues when reprojecting,
+        # cf. apps/ogr2ogr_lib.cpp:ApplySpatialFilter()
+        segment_distance_metre = 10 * 1000
+        if geom_srs.IsGeographic():
+            cloned = True
+            geom = geom.Clone()
+            dfMaxLength = segment_distance_metre / math.radians(geom_srs.GetSemiMajor())
+            geom.Segmentize(dfMaxLength)
+        elif geom_srs.IsProjected():
+            cloned = True
+            geom = geom.Clone()
+            dfMaxLength = segment_distance_metre / geom_srs.GetLinearUnits()
+            geom.Segmentize(dfMaxLength)
+
+    if geom_srs.IsSame(srs):
+        return geom
+
+    if not cloned:
+        geom = geom.Clone()
+    if geom.TransformTo(srs) != ogr.OGRERR_NONE:
+        raise RuntimeError(f'Could not transform {geom.ExportToWkt()} to {srs.GetName()}')
+    return geom
 
 def formatTZFlag(tzFlag : int) -> str:
     """Pretty-print timezone flag, cf. ogr/ogrutils.cpp:OGRGetISO8601DateTime()"""