From 074d71701d3cc423f7f3ba653fea3eb92dba269b Mon Sep 17 00:00:00 2001
From: Guilhem Moulin <guilhem@fripost.org>
Date: Mon, 10 Jun 2024 23:43:12 +0200
Subject: Add support for reprojection into the destination SRS.

The configured extent is always expressed in the destination SRS, so it
needs to be transformed into the source SRS.  Like apps/ogr2ogr_lib.cpp,
we segmentize it to make sure it is sufficiently densified.
---
 webmap-import | 73 +++++++++++++++++++++++++++++++++++++++++++++++++++--------
 1 file changed, 63 insertions(+), 10 deletions(-)

diff --git a/webmap-import b/webmap-import
index 2939ca1..39ca651 100755
--- a/webmap-import
+++ b/webmap-import
@@ -23,6 +23,7 @@ import logging
 import argparse
 import tempfile
 import re
+import math
 from fnmatch import fnmatchcase
 from pathlib import Path
 
@@ -430,6 +431,28 @@ def getSRS(srs_str):
     logging.debug('Default SRS: "%s" (%s)', srs.ExportToProj4(), srs.GetName())
     return srs
 
+def getExtent(extent):
+    if extent is None:
+        return
+
+    if not (isinstance(extent, list) or isinstance(extent, tuple)) or len(extent) != 4:
+        raise Exception(f'Invalid extent {extent}')
+
+    logging.debug('Configured extent in %s: %s',
+        srs.GetName(), ', '.join(map(str, extent)))
+
+    ring = ogr.Geometry(ogr.wkbLinearRing)
+    ring.AddPoint_2D(extent[0], extent[1])
+    ring.AddPoint_2D(extent[2], extent[1])
+    ring.AddPoint_2D(extent[2], extent[3])
+    ring.AddPoint_2D(extent[0], extent[3])
+    ring.AddPoint_2D(extent[0], extent[1])
+
+    polygon = ogr.Geometry(ogr.wkbPolygon)
+    polygon.AddGeometry(ring)
+
+    return polygon
+
 # Validate the output layer against the provided SRS and creation options
 def validateOutputLayer(lyr, srs=None, options=None):
     ok = True
@@ -842,6 +865,22 @@ def importSource2(lyr_dst, path, args={}, basedir=None, extent=None):
     logging.info('Importing layer %s from "%s"', msg, path)
     canIgnoreFields = lyr.TestCapability(ogr.OLCIgnoreFields)
 
+    srs = lyr.GetSpatialRef()
+    srs_dst = lyr_dst.GetSpatialRef()
+    if srs_dst.IsSame(srs):
+        logging.debug('Both source and destination have the same SRS (%s), skipping transformation',
+            srs_dst.GetName())
+        ct = None
+    else:
+        # TODO Use GetSupportedSRSList() and SetActiveSRS() with GDAL ≥3.7.0
+        # when possible, see apps/ogr2ogr_lib.cpp
+        logging.debug('Creating transforming from source SRS (%s) to destination SRS (%s)',
+            srs.GetName(), srs_dst.GetName())
+        ct = osr.CoordinateTransformation(srs, srs_dst)
+        if ct is None:
+            raise Exception(f'Could not create transformation from source SRS ({srs.GetName()}) '
+                + f'to destination SRS ({srs_dst.GetName()})')
+
     defn = lyr.GetLayerDefn()
     n = defn.GetGeomFieldCount()
     if n != 1: # XXX
@@ -875,8 +914,25 @@ def importSource2(lyr_dst, path, args={}, basedir=None, extent=None):
 
     count1 = -1
     if args.get('spatial-filter', True) and extent is not None:
-        # TODO reprojection
-        lyr.SetSpatialFilterRect(*extent)
+        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 extent2.TransformTo(srs) != ogr.OGRERR_NONE:
+                raise Exception(f'Could not transform extent {extent.ExportToWkt()} to {srs.GetName()}')
+
+        #logging.debug('Applying extent: %s', extent2.ExportToWkt())
+        lyr.SetSpatialFilter(extent2)
+
         if lyr.TestCapability(ogr.OLCFastFeatureCount):
             count1 = lyr.GetFeatureCount(force=0)
 
@@ -887,9 +943,6 @@ def importSource2(lyr_dst, path, args={}, basedir=None, extent=None):
         else:
             logging.info('Source layer "%s" has %d features', layername, count0)
 
-    #print(lyr.GetSupportedSRSList())
-    # TODO reprojection
-
     defn_dst = lyr_dst.GetLayerDefn()
     eGType_dst = defn_dst.GetGeomType()
 
@@ -897,12 +950,14 @@ def importSource2(lyr_dst, path, args={}, basedir=None, extent=None):
     #mismatch = 0
     feature = lyr.GetNextFeature()
     while feature is not None:
-        #print(lyr.GetFeaturesRead())
         feature2 = ogr.Feature(defn_dst)
         # TODO MakeValid, force geometry, dimension
         feature2.SetFromWithMap(feature, False, fieldMap)
 
         geom = feature2.GetGeometryRef()
+        if ct is not None and geom.Transform(ct) != ogr.OGRERR_NONE:
+            raise Exception('Could not apply coordinate transformation')
+
         eGType = geom.GetGeometryType()
         if eGType != eGType_dst:
             if eGType == ogr.wkbPolygon and eGType_dst == ogr.wkbMultiPolygon:
@@ -978,10 +1033,7 @@ if __name__ == '__main__':
 
     # get configured Spatial Reference System and extent
     srs = getSRS(common.config.get('SRS', None))
-    extent = common.config.get('extent', None)
-    if extent is not None:
-        logging.debug('Configured extent in %s: %s',
-            srs.GetName(), ', '.join(map(str, extent)))
+    extent = getExtent(common.config.get('extent', None))
 
     cachedir = Path(args.cachedir) if args.cachedir is not None else None
     rv = 0
@@ -1062,4 +1114,5 @@ if __name__ == '__main__':
 
     dso = None
     srs = None
+    extent = None
     exit(rv)
-- 
cgit v1.2.3