aboutsummaryrefslogtreecommitdiffstats
path: root/webmap-import
diff options
context:
space:
mode:
authorGuilhem Moulin <guilhem@fripost.org>2024-06-10 23:43:12 +0200
committerGuilhem Moulin <guilhem@fripost.org>2024-06-11 02:08:56 +0200
commit074d71701d3cc423f7f3ba653fea3eb92dba269b (patch)
treec06fd74702b6754364634f7e91983730c2ceeda4 /webmap-import
parent1e21ccecea0b81f13cea220b69e387542a861551 (diff)
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.
Diffstat (limited to 'webmap-import')
-rwxr-xr-xwebmap-import73
1 files 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)