aboutsummaryrefslogtreecommitdiffstats
path: root/webmap-import
diff options
context:
space:
mode:
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)