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.

2 files changed, 44 insertions, 27 deletions

diff --git a/config.yml b/config.yml
index 89a32c7..bd9dfc5 100644
--- a/config.yml
+++ b/config.yml
@@ -2,19 +2,21 @@
 # Spatial Reference System
 SRS: 'EPSG:3006' # a.k.a. SWEREF99 TM, cf. https://epsg.io/3006
 
-# (2D) extent
+# (2D) extent: minX, minY, maxX, maxY
+#
+# Lantmäteriet uses a tile-scheme where the origin (upper-left corner) is at N8500000
+# E-1200000 (SWEREF99 TM), where each tile is 256×256 pixels, and where the resolution at
+# level 0 is 4096m per pixel.
+#
+#   https://www.lantmateriet.se/globalassets/geodata/geodatatjanster/tb_twk_visning_cache_v1.1.0.pdf
+#   https://www.lantmateriet.se/globalassets/geodata/geodatatjanster/tb_twk_visning-oversiktlig_v1.0.3.pdf
+#
+# We set the extent to a 4×4 tiles square at level 2 somehow centered on Norrbotten and
+# Västerbotten.  This represent a TILEROW (x) offset of 5, and a TILECOL (y) offset of 2.
+# A 4×8 tiles rectangle with the same upper-left and upper-right coordinates can be used
+# to cover the entire country.
+#
 extent:
-  # Lantmäteriet uses a tile-scheme where the origin (upper-left corner) is at N8500000
-  # E-1200000 (SWEREF99 TM), where each tile is 256×256 pixels, and where the resolution
-  # at level 0 is 4096m per pixel.
-  #
-  #   https://www.lantmateriet.se/globalassets/geodata/geodatatjanster/tb_twk_visning_cache_v1.1.0.pdf
-  #   https://www.lantmateriet.se/globalassets/geodata/geodatatjanster/tb_twk_visning-oversiktlig_v1.0.3.pdf
-  #
-  # We set the extent to a 4×4 tiles square at level 2 somehow centered on Norrbotten and
-  # Västerbotten.  This represent a TILEROW (x) offset of 5, and a TILECOL (y) offset of 2.
-  # A 4×8 tiles rectangle with the same upper-left and upper-right coordinates can be used
-  # to cover the entire country.
   -  110720
   - 6927136 # alternatively 5878560 for the entire country
   - 1159296
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