Fix extent logic when the SRS of the output layer is not the destination SRS.

The extent is expressed in config['SRS'] in traditional GIS order
(easting/northing ordered: minX, minY, maxX, maxY), but the destination
layers might be pre-existing and use other SRS:es or mapping strategy.

1 files changed, 30 insertions, 15 deletions

diff --git a/webmap-import b/webmap-import
index 0d65db3..3242b3a 100755
--- a/webmap-import
+++ b/webmap-import
@@ -431,12 +431,18 @@ def getSRS(srs_str):
     logging.debug('Default SRS: "%s" (%s)', srs.ExportToProj4(), srs.GetName())
     return srs
 
-def getExtent(extent):
+# Convert extent [minX, minY, maxX, maxY] into a polygon and assign the
+# given SRS.  Like apps/ogr2ogr_lib.cpp, we segmentize the polygon to
+# make sure it is sufficiently densified when transforming to source
+# layer SRS for spatial filtering.
+def getExtent(extent, srs=None):
     if extent is None:
         return
 
     if not (isinstance(extent, list) or isinstance(extent, tuple)) or len(extent) != 4:
         raise Exception(f'Invalid extent {extent}')
+    elif srs is None:
+        raise Exception('Configured extent but no SRS')
 
     logging.debug('Configured extent in %s: %s',
         srs.GetName(), ', '.join(map(str, extent)))
@@ -451,6 +457,21 @@ def getExtent(extent):
     polygon = ogr.Geometry(ogr.wkbPolygon)
     polygon.AddGeometry(ring)
 
+    # we expressed extent as minX, minY, maxX, maxY (easting/northing
+    # ordered, i.e., in traditional GIS order)
+    srs2 = srs.Clone()
+    srs2.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    polygon.AssignSpatialReference(srs2)
+    polygon.TransformTo(srs)
+
+    segment_distance_metre = 10 * 1000
+    if srs.IsGeographic():
+        dfMaxLength = segment_distance_metre / math.radians(srs.GetSemiMajor())
+        polygon.Segmentize(dfMaxLength)
+    elif srs.IsProjected():
+        dfMaxLength = segment_distance_metre / srs.GetLinearUnits()
+        polygon.Segmentize(dfMaxLength)
+
     return polygon
 
 # Validate the output layer against the provided SRS and creation options
@@ -866,6 +887,9 @@ def importSource2(lyr_dst, path, args={}, basedir=None, extent=None):
     canIgnoreFields = lyr.TestCapability(ogr.OLCIgnoreFields)
 
     srs = lyr.GetSpatialRef()
+    if srs is None:
+        raise Exception('Source layer has no SRS')
+
     srs_dst = lyr_dst.GetSpatialRef()
     if srs_dst.IsSame(srs):
         logging.debug('Both source and destination have the same SRS (%s), skipping coordinate transformation',
@@ -914,19 +938,10 @@ def importSource2(lyr_dst, path, args={}, basedir=None, extent=None):
 
     count1 = -1
     if args.get('spatial-filter', True) and extent is not None:
-        extent2 = extent.Clone()
-        extent2.AssignSpatialReference(srs_dst)
-
-        if srs is not None and not srs.IsSame(srs_dst):
-            # make sure to densify it sufficiently to avoid issues, cf. apps/ogr2ogr_lib.cpp
-            segment_distance_metre = 10 * 1000
-            if srs_dst.IsGeographic():
-                dfMaxLength = segment_distance_metre / math.radians(srs_dst.GetSemiMajor())
-                extent2.Segmentize(dfMaxLength)
-            elif srs_dst.IsProjected():
-                dfMaxLength = segment_distance_metre / srs_dst.GetLinearUnits()
-                extent2.Segmentize(dfMaxLength)
-
+        if extent.GetSpatialReference().IsSame(srs):
+            extent2 = extent
+        else:
+            extent2 = extent.Clone()
             if extent2.TransformTo(srs) != ogr.OGRERR_NONE:
                 raise Exception(f'Could not transform extent {extent.ExportToWkt()} to {srs.GetName()}')
 
@@ -1033,7 +1048,7 @@ if __name__ == '__main__':
 
     # get configured Spatial Reference System and extent
     srs = getSRS(common.config.get('SRS', None))
-    extent = getExtent(common.config.get('extent', None))
+    extent = getExtent(common.config.get('extent', None), srs=srs)
 
     cachedir = Path(args.cachedir) if args.cachedir is not None else None
     rv = 0