diff options
| -rwxr-xr-x | webmap-import | 73 |
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) |