diff options
| -rw-r--r-- | common_gdal.py | 60 | ||||
| -rw-r--r-- | import_source.py | 21 | ||||
| -rwxr-xr-x | webmap-import | 4 |
3 files changed, 53 insertions, 32 deletions
diff --git a/common_gdal.py b/common_gdal.py index 4fea54a..35511a3 100644 --- a/common_gdal.py +++ b/common_gdal.py @@ -94,15 +94,11 @@ def getSRS(srs_str : Optional[str]) -> osr.SpatialReference: return srs def getExtent(extent : Optional[tuple[float, float, float, float]], - srs : Optional[osr.SpatialReference] = None) -> tuple[ogr.Geometry, ogr.Geometry]: + srs : Optional[osr.SpatialReference] = None) -> ogr.Geometry: """Convert extent (minX, minY, maxX, maxY) into a polygon and assign the - given SRS. Return a pair with the densified and non-densified extent. - Like apps/ogr2ogr_lib.cpp, the former is obtained by segmentizing the - polygon to make sure it is sufficiently densified when transforming to - source layer SRS for spatial filtering.""" - + given SRS.""" if extent is None: - return None, None + return None if not isinstance(extent, tuple) or len(extent) != 4: raise RuntimeError(f'Invalid extent {extent}') @@ -130,18 +126,44 @@ def getExtent(extent : Optional[tuple[float, float, float, float]], if not srs2.IsSame(srs): polygon.TransformTo(srs) - # densify the rectangle to avoid issues when reprojecting to the - # source layer SRS, cf. apps/ogr2ogr_lib.cpp:ApplySpatialFilter() - polygon_dense = polygon.Clone() - segment_distance_metre = 10 * 1000 - if srs.IsGeographic(): - dfMaxLength = segment_distance_metre / math.radians(srs.GetSemiMajor()) - polygon_dense.Segmentize(dfMaxLength) - elif srs.IsProjected(): - dfMaxLength = segment_distance_metre / srs.GetLinearUnits() - polygon_dense.Segmentize(dfMaxLength) - - return polygon_dense, polygon + return polygon + +def getSpatialFilterFromGeometry(geom : ogr.Geometry, + srs : osr.SpatialReference, + fail_if_srs_not_equivalent : bool = False) -> ogr.Geometry: + """Make the geometry suitable for use as a spatial filter. It is + densified the SRS:s are not equivalent, unless fail_if_srs_not_equivalent is True in + which case it errors out (this ensure clipping geometries remain simple).""" + cloned = False + geom_srs = geom.GetSpatialReference() + if not geom_srs.IsSame(srs, [ 'IGNORE_DATA_AXIS_TO_SRS_AXIS_MAPPING=YES', + 'CRITERION=EQUIVALENT']): + if fail_if_srs_not_equivalent: + raise RuntimeError(f'Geometry {geom.ExportToWkt()} has SRS ' + + geom_srs.ExportToPrettyWkt() + ',\nexpected ' + + srs.ExportToPrettyWkt()) + # densify the geometry (a rectangle) to avoid issues when reprojecting, + # cf. apps/ogr2ogr_lib.cpp:ApplySpatialFilter() + segment_distance_metre = 10 * 1000 + if geom_srs.IsGeographic(): + cloned = True + geom = geom.Clone() + dfMaxLength = segment_distance_metre / math.radians(geom_srs.GetSemiMajor()) + geom.Segmentize(dfMaxLength) + elif geom_srs.IsProjected(): + cloned = True + geom = geom.Clone() + dfMaxLength = segment_distance_metre / geom_srs.GetLinearUnits() + geom.Segmentize(dfMaxLength) + + if geom_srs.IsSame(srs): + return geom + + if not cloned: + geom = geom.Clone() + if geom.TransformTo(srs) != ogr.OGRERR_NONE: + raise RuntimeError(f'Could not transform {geom.ExportToWkt()} to {srs.GetName()}') + return geom def formatTZFlag(tzFlag : int) -> str: """Pretty-print timezone flag, cf. ogr/ogrutils.cpp:OGRGetISO8601DateTime()""" diff --git a/import_source.py b/import_source.py index 6fe4acf..04dea4b 100644 --- a/import_source.py +++ b/import_source.py @@ -43,7 +43,13 @@ from osgeo.gdalconst import ( from osgeo import gdalconst from common import BadConfiguration, escape_identifier, escape_literal_str -from common_gdal import gdalSetOpenExArgs, gdalGetMetadataItem, gdalVersionMin, formatTZFlag +from common_gdal import ( + gdalSetOpenExArgs, + gdalGetMetadataItem, + gdalVersionMin, + formatTZFlag, + getSpatialFilterFromGeometry, +) def openOutputDS(def_dict : dict[str, Any]) -> gdal.Dataset: """Open and return the output DS. It is created if create=False or @@ -692,16 +698,9 @@ def _importSource2(lyr_dst : ogr.Layer, path : str, args : dict[str,Any], count1 = -1 if args.get('spatial-filter', True) and extent is not None: - if extent.GetSpatialReference().IsSame(srs): - extent2 = extent - else: - extent2 = extent.Clone() - if extent2.TransformTo(srs) != ogr.OGRERR_NONE: - raise RuntimeError(f'Could not transform extent {extent.ExportToWkt()} ' - f'to {srs.GetName()}') - - #logging.debug('Applying extent: %s', extent2.ExportToWkt()) - lyr.SetSpatialFilter(extent2) + spatialFilter = getSpatialFilterFromGeometry(extent, srs) + logging.debug('Setting spatial filter to %s', spatialFilter.ExportToWkt()) + lyr.SetSpatialFilter(spatialFilter) if lyr.TestCapability(ogr.OLCFastFeatureCount): count1 = lyr.GetFeatureCount(force=0) diff --git a/webmap-import b/webmap-import index f20fdef..dafff78 100755 --- a/webmap-import +++ b/webmap-import @@ -613,7 +613,7 @@ def main() -> NoReturn: # get configured Spatial Reference System and extent srs = getSRS(config.get('SRS', None)) - extent = getExtent(config.get('extent', None), srs=srs)[0] + extent = getExtent(config.get('extent', None), srs=srs) if args.lockfile is not None: # obtain an exclusive lock and don't release it until exiting the program @@ -711,8 +711,8 @@ def main() -> NoReturn: finally: lyr_cache = None dso = None - srs = None extent = None + srs = None sys.exit(rv) gdal.UseExceptions() |