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.

3 files changed, 53 insertions, 32 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()"""
diff --git a/import_source.py b/import_source.py
index 6fe4acf..04dea4b 100644
--- a/import_source.py
+++ b/import_source.py
@@ -43,7 +43,13 @@ from osgeo.gdalconst import (
 from osgeo import gdalconst
 
 from common import BadConfiguration, escape_identifier, escape_literal_str
-from common_gdal import gdalSetOpenExArgs, gdalGetMetadataItem, gdalVersionMin, formatTZFlag
+from common_gdal import (
+    gdalSetOpenExArgs,
+    gdalGetMetadataItem,
+    gdalVersionMin,
+    formatTZFlag,
+    getSpatialFilterFromGeometry,
+)
 
 def openOutputDS(def_dict : dict[str, Any]) -> gdal.Dataset:
     """Open and return the output DS.  It is created if create=False or
@@ -692,16 +698,9 @@ def _importSource2(lyr_dst : ogr.Layer, path : str, args : dict[str,Any],
 
     count1 = -1
     if args.get('spatial-filter', True) and extent is not None:
-        if extent.GetSpatialReference().IsSame(srs):
-            extent2 = extent
-        else:
-            extent2 = extent.Clone()
-            if extent2.TransformTo(srs) != ogr.OGRERR_NONE:
-                raise RuntimeError(f'Could not transform extent {extent.ExportToWkt()} '
-                                   f'to {srs.GetName()}')
-
-        #logging.debug('Applying extent: %s', extent2.ExportToWkt())
-        lyr.SetSpatialFilter(extent2)
+        spatialFilter = getSpatialFilterFromGeometry(extent, srs)
+        logging.debug('Setting spatial filter to %s', spatialFilter.ExportToWkt())
+        lyr.SetSpatialFilter(spatialFilter)
 
         if lyr.TestCapability(ogr.OLCFastFeatureCount):
             count1 = lyr.GetFeatureCount(force=0)
diff --git a/webmap-import b/webmap-import
index f20fdef..dafff78 100755
--- a/webmap-import
+++ b/webmap-import
@@ -613,7 +613,7 @@ def main() -> NoReturn:
 
     # get configured Spatial Reference System and extent
     srs = getSRS(config.get('SRS', None))
-    extent = getExtent(config.get('extent', None), srs=srs)[0]
+    extent = getExtent(config.get('extent', None), srs=srs)
 
     if args.lockfile is not None:
         # obtain an exclusive lock and don't release it until exiting the program
@@ -711,8 +711,8 @@ def main() -> NoReturn:
     finally:
         lyr_cache = None
         dso = None
-        srs = None
         extent = None
+        srs = None
     sys.exit(rv)
 
 gdal.UseExceptions()