diff options
| -rw-r--r-- | common_gdal.py | 11 | ||||
| -rw-r--r-- | export_mvt.py | 29 |
2 files changed, 7 insertions, 33 deletions
diff --git a/common_gdal.py b/common_gdal.py index 35511a3..7333f58 100644 --- a/common_gdal.py +++ b/common_gdal.py @@ -128,20 +128,13 @@ def getExtent(extent : Optional[tuple[float, float, float, float]], return polygon -def getSpatialFilterFromGeometry(geom : ogr.Geometry, - srs : osr.SpatialReference, - fail_if_srs_not_equivalent : bool = False) -> ogr.Geometry: +def getSpatialFilterFromGeometry(geom : ogr.Geometry, srs : osr.SpatialReference) -> 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).""" + densified the SRS:s are not equivalent.""" 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 diff --git a/export_mvt.py b/export_mvt.py index e9b54f4..73d2ac7 100644 --- a/export_mvt.py +++ b/export_mvt.py @@ -164,12 +164,9 @@ def exportSourceLayer(ds_src : gdal.Dataset, if extent is None: spatialFilter = None - clip = None else: # transform extent to source SRS spatialFilter = getSpatialFilterFromGeometry(extent, srs_src) - # ensure clipping geometries remain simple (mere rectangles) - clip = getSpatialFilterFromGeometry(extent, srs_dst, fail_if_srs_not_equivalent = True) transform_geometry = layerdef.get('transform-geometry', None) columns = [ 'm.' + escape_identifier(lyr_src.GetFIDColumn()) ] @@ -198,45 +195,29 @@ def exportSourceLayer(ds_src : gdal.Dataset, logging.debug('Source layer "%s" has %d features, of which %d are to be exported', layername, count0, count1) + feature_count = 0 defn_dst = lyr_dst.GetLayerDefn() - eGType_dst = defn_dst.GetGeomType() - dGeomIsUnknown = ogr.GT_Flatten(eGType_dst) == ogr.wkbUnknown - - count = count_degenerates = 0 feature = lyr_src.GetNextFeature() while feature is not None: geom = feature.GetGeometryRef().Clone() if ct is not None and geom.Transform(ct) != ogr.OGRERR_NONE: raise RuntimeError('Could not apply coordinate transformation') geom.FlattenTo2D() - if clip is not None: - # GDAL generates tiles beyond the tiling scheme extents so we need to - # clip manually - dim = geom.GetDimension() - geom = geom.Intersection(clip) - if geom.GetDimension() != dim: - # degenerate intersection, skip feature - count_degenerates += 1 - feature = lyr_src.GetNextFeature() - continue - if not dGeomIsUnknown: - geom = ogr.ForceTo(geom, eGType_dst) - feature2 = ogr.Feature(defn_dst) feature2.SetGeometryDirectly(geom) feature2.SetFID(feature.GetFID()) if lyr_dst.CreateFeature(feature2) != ogr.OGRERR_NONE: raise RuntimeError(f'Could not transfer source feature #{feature.GetFID()}') - count += 1 + feature_count += 1 feature = lyr_src.GetNextFeature() - logging.info('Exported %d features (%d degenerate skipped) to MVT layer "%s" from "%s"', - count, count_degenerates, lyr_dst.GetName(), layername) + logging.info('Exported %d features to MVT layer "%s" from "%s"', + feature_count, lyr_dst.GetName(), layername) finally: ds_src.ReleaseResultSet(lyr_src) lyr_src = None - return count + return feature_count def list_files(top : str, dir_fd : Optional[int] = None, |