MVT: Don't clip features manually.

It's done automatically by the MVT driver.  Also we don't want to clip
at the extent boundaries, but instead leave som headroom so the clipped
polygon border is not visible on the map.  The MVT driver does that.

It seems that GDAL 3.6.2 from Debian Bookworm generates too many tiles
though.  It yields the following tile counts for group ‘ren’:

    no manual clipping, BUFFER=32: 83718 tiles  [min=33 B, max=117.70 kiB, sum=15.73 MiB, avg=197 B]
    no manual clipping, BUFFER=0:  83676 tiles
    clip at extent,     BUFFER=32: 76256 tiles

GDAL 3.10.3 from Debian Trixie yields less surprising tile counts:

    no manual clipping, BUFFER=32: 75972 tiles [min=33 B, max=128.16 kiB, sum=15.10 MiB, avg=208 B]
    no manual clipping, BUFFER=0:  75939 tiles
    clip at extent,     BUFFER=32: 75972 tiles

(Interesting to see that the largest tile — 0/0/0.pbf — is over 10kiB
larger with the more recent GDAL version, also.)

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,