Fix column/polygon_points sync on degenerate polygon drop (#1151)

_simplify_polygons now accepts and returns column alongside
polygon_points so both stay in sync when degenerate polygons
are dropped. Added tests for holes, transform, single-pixel
raster, and column/pp length consistency.

Author: brendancol

xrspatial/polygonize.py

@@ -1272,7 +1272,8 @@ def _chain_key(chain):
     return (start, end, interior)
 
 
-def _simplify_polygons(polygon_points, tolerance, method="douglas-peucker"):
+def _simplify_polygons(column, polygon_points, tolerance,
+                       method="douglas-peucker"):
     """Topology-preserving simplification of all polygons.
 
     Uses shared-edge decomposition: finds junction vertices, splits
@@ -1281,6 +1282,8 @@ def _simplify_polygons(polygon_points, tolerance, method="douglas-peucker"):
 
     Parameters
     ----------
+    column : list
+        Pixel values corresponding to each polygon.
     polygon_points : list of list of np.ndarray
         Output of polygonize backend: list of polygons, each polygon
         is [exterior_ring, *hole_rings].
@@ -1292,11 +1295,12 @@ def _simplify_polygons(polygon_points, tolerance, method="douglas-peucker"):
 
     Returns
     -------
-    list of list of np.ndarray
-        Same structure as input, with simplified coordinates.
+    (list, list of list of np.ndarray)
+        Filtered column and simplified polygon_points.  Polygons whose
+        exterior ring collapses below 4 vertices are dropped from both.
     """
     if tolerance <= 0:
-        return polygon_points
+        return column, polygon_points
 
     # Step 1: Find junctions.
     junctions = _find_junctions(polygon_points)
@@ -1330,6 +1334,7 @@ def _simplify_polygons(polygon_points, tolerance, method="douglas-peucker"):
         ring_info.append(poly_info)
 
     # Step 4: Reassemble rings.
+    result_column = []
     result = []
     for poly_idx, rings in enumerate(polygon_points):
         new_rings = []
@@ -1371,9 +1376,10 @@ def _simplify_polygons(polygon_points, tolerance, method="douglas-peucker"):
         for ring in new_rings[1:]:
             if len(ring) >= 4:
                 filtered.append(ring)
+        result_column.append(column[poly_idx])
         result.append(filtered)
 
-    return result
+    return result_column, result
 
 
 def _merge_polygon_rings(polys_list):
@@ -1608,8 +1614,9 @@ def polygonize(
 
     # Apply simplification if requested.
     if simplify_tolerance is not None and simplify_tolerance > 0:
-        polygon_points = _simplify_polygons(
-            polygon_points, simplify_tolerance, method=simplify_method)
+        column, polygon_points = _simplify_polygons(
+            column, polygon_points, simplify_tolerance,
+            method=simplify_method)
 
     # Convert to requested return_type.
     if return_type == "numpy":

xrspatial/tests/test_polygonize.py

@@ -593,7 +593,8 @@ def test_simplify_polygons_reduces_vertices(self):
         data = xr.DataArray(raster)
         column_orig, pp_orig = polygonize(data, return_type="numpy")
 
-        pp_simplified = _simplify_polygons(pp_orig, tolerance=1.5)
+        _, pp_simplified = _simplify_polygons(
+            column_orig, pp_orig, tolerance=1.5)
 
         # Vertex count should be reduced for at least one polygon.
         orig_total_verts = sum(
@@ -624,7 +625,7 @@ def test_simplify_polygons_topology_preserved(self):
         data = xr.DataArray(raster)
         column, pp = polygonize(data, return_type="numpy")
 
-        pp_simplified = _simplify_polygons(pp, tolerance=0.0)
+        _, pp_simplified = _simplify_polygons(column, pp, tolerance=0.0)
 
         # With tolerance=0, no vertices should be removed; output
         # should match input exactly.
@@ -648,7 +649,7 @@ def test_simplify_polygons_shared_edge_identical(self):
         data = xr.DataArray(raster)
         column, pp = polygonize(data, return_type="numpy")
 
-        pp_simplified = _simplify_polygons(pp, tolerance=1.5)
+        _, pp_simplified = _simplify_polygons(column, pp, tolerance=1.5)
 
         # Check total area is preserved (which requires no gaps/overlaps).
         total_orig = sum(
@@ -708,12 +709,28 @@ def test_simplify_polygons_drops_degenerate(self):
 
         # Find the small polygon (value 1, single pixel).
         small_idx = column.index(1)
-        small_poly = [pp[small_idx]]  # wrap in list for _simplify_polygons
+        small_col = [column[small_idx]]
+        small_poly = [pp[small_idx]]
 
         # Large tolerance should collapse the single-pixel polygon.
-        result = _simplify_polygons(small_poly, tolerance=10.0)
+        result_col, result_pp = _simplify_polygons(
+            small_col, small_poly, tolerance=10.0)
+        # Column and polygon_points must stay in sync.
+        assert len(result_col) == len(result_pp)
         # The polygon should be dropped since its exterior collapses.
-        assert len(result) == 0 or all(len(r[0]) >= 4 for r in result)
+        assert len(result_pp) == 0 or all(len(r[0]) >= 4 for r in result_pp)
+
+    def test_simplify_column_pp_sync_through_public_api(self):
+        """Public API must keep column and polygon_points in sync."""
+        # Raster with a single-pixel region that may collapse.
+        raster = np.array([
+            [1, 2, 2],
+            [2, 2, 2],
+            [2, 2, 2],
+        ], dtype=np.int64)
+        data = xr.DataArray(raster)
+        col, pp = polygonize(data, simplify_tolerance=10.0)
+        assert len(col) == len(pp)
 
 
 class TestPolygonizeSimplify:
@@ -879,6 +896,45 @@ def test_simplify_with_awkward(self):
                             return_type="awkward")
         assert result is not None
 
+    def test_simplify_with_holes(self):
+        """Simplification should handle polygons with interior holes."""
+        raster = np.zeros((8, 8), dtype=np.int32)
+        raster[1:7, 1:7] = 1
+        raster[3:5, 3:5] = 0  # hole in the 1-region
+        data = xr.DataArray(raster)
+
+        col_orig, pp_orig = polygonize(data)
+        col_simp, pp_simp = polygonize(data, simplify_tolerance=0.5)
+
+        assert len(col_simp) == len(pp_simp)
+        # Area should be preserved.
+        total_orig = sum(
+            sum(calc_boundary_area(r) for r in rings) for rings in pp_orig)
+        total_simp = sum(
+            sum(calc_boundary_area(r) for r in rings) for rings in pp_simp)
+        assert_allclose(total_simp, total_orig, atol=1e-10)
+
+    def test_simplify_with_transform(self):
+        """Simplification should work with affine transform."""
+        raster = np.array([
+            [1, 1, 2, 2],
+            [1, 1, 2, 2],
+            [1, 1, 2, 2],
+        ], dtype=np.int64)
+        transform = np.array([10.0, 0.0, 100.0, 0.0, 10.0, 200.0])
+        data = xr.DataArray(raster)
+        col, pp = polygonize(data, transform=transform,
+                             simplify_tolerance=5.0)
+        assert len(col) == len(pp)
+        assert len(col) == 2
+
+    def test_simplify_single_pixel_raster(self):
+        """Simplification on a 1x1 raster should not crash."""
+        raster = np.array([[5]], dtype=np.int64)
+        data = xr.DataArray(raster)
+        col, pp = polygonize(data, simplify_tolerance=1.0)
+        assert len(col) == len(pp)
+
 
 @pytest.mark.skipif(da is None, reason="dask not installed")
 class TestPolygonizeSimplifyDask: