gdalwarp: cutline zero-width sliver enhancement: avoid producing inva…

…lid polygons
rouault · Feb 18, 2024 · 50bfe4d · 50bfe4d
1 parent 7500c17
commit 50bfe4d
Show file tree

Hide file tree

Showing 2 changed files with 129 additions and 4 deletions.
diff --git a/apps/gdalwarp_lib.cpp b/apps/gdalwarp_lib.cpp
@@ -4204,16 +4204,56 @@ static void RemoveZeroWidthSlivers(OGRGeometry *poGeom)
     const OGRwkbGeometryType eType = wkbFlatten(poGeom->getGeometryType());
     if (eType == wkbMultiPolygon)
     {
-        for (auto poSubGeom : *(poGeom->toMultiPolygon()))
+        auto poMP = poGeom->toMultiPolygon();
+        int nNumGeometries = poMP->getNumGeometries();
+        for (int i = 0; i < nNumGeometries; /* incremented in loop */)
         {
-            RemoveZeroWidthSlivers(poSubGeom);
+            auto poPoly = poMP->getGeometryRef(i);
+            RemoveZeroWidthSlivers(poPoly);
+            if (poPoly->IsEmpty())
+            {
+                CPLDebug("WARP",
+                         "RemoveZeroWidthSlivers: removing empty polygon");
+                poMP->removeGeometry(i, /* bDelete = */ true);
+                --nNumGeometries;
+            }
+            else
+            {
+                ++i;
+            }
         }
     }
     else if (eType == wkbPolygon)
     {
-        for (auto poSubGeom : *(poGeom->toPolygon()))
+        auto poPoly = poGeom->toPolygon();
+        if (auto poExteriorRing = poPoly->getExteriorRing())
         {
-            RemoveZeroWidthSlivers(poSubGeom);
+            RemoveZeroWidthSlivers(poExteriorRing);
+            if (poExteriorRing->getNumPoints() < 4)
+            {
+                poPoly->empty();
+                return;
+            }
+        }
+        int nNumInteriorRings = poPoly->getNumInteriorRings();
+        for (int i = 0; i < nNumInteriorRings; /* incremented in loop */)
+        {
+            auto poRing = poPoly->getInteriorRing(i);
+            RemoveZeroWidthSlivers(poRing);
+            if (poRing->getNumPoints() < 4)
+            {
+                CPLDebug(
+                    "WARP",
+                    "RemoveZeroWidthSlivers: removing empty interior ring");
+                constexpr int OFFSET_EXTERIOR_RING = 1;
+                poPoly->removeRing(i + OFFSET_EXTERIOR_RING,
+                                   /* bDelete = */ true);
+                --nNumInteriorRings;
+            }
+            else
+            {
+                ++i;
+            }
         }
     }
     else if (eType == wkbLineString)

diff --git a/autotest/utilities/test_gdalwarp_lib.py b/autotest/utilities/test_gdalwarp_lib.py
@@ -30,6 +30,7 @@
 # DEALINGS IN THE SOFTWARE.
 ###############################################################################
 
+import json
 import os
 import shutil
 import struct
@@ -2595,6 +2596,90 @@ def test_gdalwarp_lib_cutline_zero_width_sliver():
     assert ds is not None
 
 
+###############################################################################
+# Test cutline with zero-width sliver
+
+
+def test_gdalwarp_lib_cutline_zero_width_sliver_remove_empty_polygon():
+
+    geojson = {
+        "type": "MultiPolygon",
+        "coordinates": [
+            [
+                [
+                    [-101.43346, 36.91886],
+                    [-101.43337, 36.91864],
+                    [-101.43332, 36.91865],
+                    [-101.43342, 36.91888],
+                    [-101.43346, 36.91886],
+                ]
+            ],
+            # The below polygon has a zero-width sliver
+            [
+                [
+                    [-101.4311, 36.91909],
+                    [-101.43106, 36.91913],
+                    [-101.43111, 36.91908],
+                    [-101.4311, 36.91909],
+                ]
+            ],
+        ],
+    }
+    gdal.FileFromMemBuffer("/vsimem/cutline.geojson", json.dumps(geojson))
+    src_ds = gdal.GetDriverByName("MEM").Create("", 100, 100)
+    srs = osr.SpatialReference()
+    srs.ImportFromEPSG(4326)
+    src_ds.SetSpatialRef(srs)
+    src_ds.SetGeoTransform([-101.5, 0.1, 0, 37, 0, -0.1])
+    ds = gdal.Warp("", src_ds, format="MEM", cutlineDSName="/vsimem/cutline.geojson")
+    assert ds is not None
+
+
+###############################################################################
+# Test cutline with zero-width sliver
+
+
+def test_gdalwarp_lib_cutline_zero_width_sliver_remove_empty_inner_ring():
+
+    geojson = {
+        "type": "MultiPolygon",
+        "coordinates": [
+            [
+                [
+                    [-101.5, 37],
+                    [-101.4, 37],
+                    [-101.4, 36.9],
+                    [-101.5, 36.9],
+                    [-101.5, 37],
+                ],
+                # The below ring has a zero-width sliver
+                [
+                    [-101.4311, 36.91909],
+                    [-101.43106, 36.91913],
+                    [-101.43111, 36.91908],
+                    [-101.4311, 36.91909],
+                ],
+                # This one is OK
+                [
+                    [-101.49, 36.95],
+                    [-101.48, 36.95],
+                    [-101.48, 36.94],
+                    [-101.49, 36.94],
+                    [-101.49, 36.95],
+                ],
+            ]
+        ],
+    }
+    gdal.FileFromMemBuffer("/vsimem/cutline.geojson", json.dumps(geojson))
+    src_ds = gdal.GetDriverByName("MEM").Create("", 100, 100)
+    srs = osr.SpatialReference()
+    srs.ImportFromEPSG(4326)
+    src_ds.SetSpatialRef(srs)
+    src_ds.SetGeoTransform([-101.6, 0.1, 0, 37.1, 0, -0.1])
+    ds = gdal.Warp("", src_ds, format="MEM", cutlineDSName="/vsimem/cutline.geojson")
+    assert ds is not None
+
+
 ###############################################################################
 # Test support for propagating coordinate epoch