Factor out densification logic from getExtent() into own function.

And only densify if needs be.  Most sources are already in SWEREF 99
(modulo axis mapping strategy) so in pratice we can use mere rectangles
as spatial filters.

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