From a8833e35eba6a0617723ec0e751265ebc6c4bccd Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Sun, 15 Mar 2026 14:38:01 -0700
Subject: [PATCH 1/3] Improve polygonize performance: single-pass tracing, JIT
 merge helpers, batch shapely (#1008)

Replace the two-pass _follow with a single-pass implementation using a
dynamically-grown buffer. This eliminates retracing every polygon
boundary a second time, which was the dominant cost for rasters with
many small regions.

Add @ngjit to _point_in_ring, _simplify_ring, and _signed_ring_area
so the dask chunk-merge path runs compiled instead of interpreted.

Use shapely.polygons() batch constructor for hole-free polygons in
_to_geopandas (shapely 2.0+, with fallback for older versions).
---
 xrspatial/polygonize.py | 245 +++++++++++++++++++++++-----------------
 1 file changed, 141 insertions(+), 104 deletions(-)

diff --git a/xrspatial/polygonize.py b/xrspatial/polygonize.py
index 40ca8c6a..439db993 100644
--- a/xrspatial/polygonize.py
+++ b/xrspatial/polygonize.py
@@ -107,10 +107,10 @@ def _min_and_max(value0, value1):
 # facing W, and if hole is False the start is on the S edge of the pixel
 # facing E.
 #
-# There are two passes.  First pass determines the number of points and
-# allocates a numpy array big enough to take the points, the second pass saves
-# the points.  This is better than a single pass that repeatedly reallocates
-# space for the increasing number of points.
+# Single pass with a dynamically-grown buffer (doubling strategy).  Buffer
+# starts at 64 points and grows as needed.  Amortized cost is O(1) per point
+# from O(log n) reallocations, which is cheaper than retracing the boundary
+# a second time.
 #
 # Returns the region ID and the 2D array of boundary points.  The last
 # boundary point is the same as the first.
@@ -125,98 +125,104 @@ def _follow(
 ) -> Tuple[int, np.ndarray]:
     region = regions[ij]
     n = nx*ny
-    points = None  # Declare before loop for numba
 
-    for pass_ in range(2):
-        prev_forward = 0  # Start with an invalid direction.
-        if hole:
-            forward = -1  # Facing W along N edge.
-            left = -nx
-        else:
-            forward = 1  # Facing E along S edge.
-            left = nx
-
-        start_forward = forward
-        start_ij = ij
-        npoints = 0
-
-        while True:
-            if pass_ == 1:
-                if forward == 1 and not hole:
-                    # Mark pixel so that it will not be considered a future
-                    # non-hole starting pixel.
-                    visited[ij] |= 1
-                elif forward == -1 and ij+nx < n:
-                    # Mark pixel so that is will not be considered a future
-                    # hole starting pixel.
-                    visited[ij+nx] |= 2
-
-            if prev_forward != forward:
-                if pass_ == 1:
-                    # Add point.
-                    i = ij % nx
-                    j = ij // nx
-                    if forward == -1:
-                        i += 1
-                        j += 1
-                    elif forward == nx:
-                        i += 1
-                    elif forward == -nx:
-                        j += 1
-                    points[2*npoints] = i
-                    points[2*npoints+1] = j
-                npoints += 1
+    # Dynamic buffer for point coordinates (x, y pairs stored flat).
+    buf_cap = 64
+    points = np.empty(2 * buf_cap)
+    npoints = 0
+
+    prev_forward = 0  # Start with an invalid direction.
+    if hole:
+        forward = -1  # Facing W along N edge.
+        left = -nx
+    else:
+        forward = 1  # Facing E along S edge.
+        left = nx
+
+    start_forward = forward
+    start_ij = ij
+
+    while True:
+        if forward == 1 and not hole:
+            # Mark pixel so that it will not be considered a future
+            # non-hole starting pixel.
+            visited[ij] |= 1
+        elif forward == -1 and ij+nx < n:
+            # Mark pixel so that it will not be considered a future
+            # hole starting pixel.
+            visited[ij+nx] |= 2
+
+        if prev_forward != forward:
+            # Grow buffer if needed.
+            if npoints >= buf_cap:
+                old_points = points
+                buf_cap *= 2
+                points = np.empty(2 * buf_cap)
+                points[:2*npoints] = old_points[:2*npoints]
+
+            # Add point.
+            i = ij % nx
+            j = ij // nx
+            if forward == -1:
+                i += 1
+                j += 1
+            elif forward == nx:
+                i += 1
+            elif forward == -nx:
+                j += 1
+            points[2*npoints] = i
+            points[2*npoints+1] = j
+            npoints += 1
+
+        prev_forward = forward
+        ijnext = ij + forward
+        ijnext_right = ijnext - left
+
+        # Determine direction of turn.
+        if abs(forward) == 1:  # Facing E (forward == 1) or W (forward -1)
+            if _diff_row(ij, ijnext, nx):
+                turn = Turn.Left
+            elif (not _outside_domain(ijnext_right, n) and
+                    regions[ijnext_right] == region):
+                turn = Turn.Right
+            elif regions[ijnext] == region:
+                turn = Turn.Straight
+            else:
+                turn = Turn.Left
+        else:  # Facing N (forward == nx) or S (forward == -nx)
+            if _outside_domain(ijnext, n):
+                turn = Turn.Left
+            elif (not _diff_row(ijnext, ijnext_right, nx) and
+                    regions[ijnext_right] == region):
+                turn = Turn.Right
+            elif regions[ijnext] == region:
+                turn = Turn.Straight
+            else:
+                turn = Turn.Left
 
+        # Apply turn.
+        if turn == Turn.Straight:
+            ij = ijnext
+        elif turn == Turn.Left:
             prev_forward = forward
-            ijnext = ij + forward
-            ijnext_right = ijnext - left
-
-            # Determine direction of turn.
-            if abs(forward) == 1:  # Facing E (forward == 1) or W (forward -1)
-                if _diff_row(ij, ijnext, nx):
-                    turn = Turn.Left
-                elif (not _outside_domain(ijnext_right, n) and
-                        regions[ijnext_right] == region):
-                    turn = Turn.Right
-                elif regions[ijnext] == region:
-                    turn = Turn.Straight
-                else:
-                    turn = Turn.Left
-            else:  # Facing N (forward == nx) or S (forward == -nx)
-                if _outside_domain(ijnext, n):
-                    turn = Turn.Left
-                elif (not _diff_row(ijnext, ijnext_right, nx) and
-                        regions[ijnext_right] == region):
-                    turn = Turn.Right
-                elif regions[ijnext] == region:
-                    turn = Turn.Straight
-                else:
-                    turn = Turn.Left
-
-            # Apply turn.
-            if turn == Turn.Straight:
-                ij = ijnext
-            elif turn == Turn.Left:
-                prev_forward = forward
-                forward = left
-                left = -prev_forward
-            else:  # Turn.Right
-                prev_forward = forward
-                forward = -left
-                left = prev_forward
-                ij = ijnext_right
-
-            # Finished boundary when returned to start.
-            if ij == start_ij and forward == start_forward:
-                break
+            forward = left
+            left = -prev_forward
+        else:  # Turn.Right
+            prev_forward = forward
+            forward = -left
+            left = prev_forward
+            ij = ijnext_right
 
-        if pass_ == 0:
-            # End of first pass.
-            points = np.empty(2*(npoints+1))  # Note extra point at end.
+        # Finished boundary when returned to start.
+        if ij == start_ij and forward == start_forward:
+            break
 
-    points = points.reshape((-1, 2))
-    points[-1] = points[0]  # End point the same as start point.
-    return region, points
+    # Trim buffer to actual size and add closing point.
+    result = np.empty(2*(npoints+1))
+    result[:2*npoints] = points[:2*npoints]
+    result = result.reshape((-1, 2))
+    result[-1] = result[0]  # End point the same as start point.
+    return region, result
 
 
 # Generator of numba-compatible comparison functions for values.
@@ -474,11 +480,34 @@ def _to_geopandas(
     column_name: str,
 ):
     import geopandas as gpd
+    import shapely
     from shapely.geometry import Polygon
 
-    # Convert list of point arrays to shapely Polygons.
-    polygons = list(map(
-        lambda points: Polygon(points[0], points[1:]), polygon_points))
+    if hasattr(shapely, 'polygons'):
+        # Shapely 2.0+: batch-construct hole-free polygons via
+        # linearrings -> polygons pipeline (both are C-level batch ops).
+        no_holes = [i for i, pts in enumerate(polygon_points)
+                    if len(pts) == 1]
+
+        if len(no_holes) == len(polygon_points):
+            # All hole-free: batch create LinearRings then Polygons.
+            rings = [shapely.linearrings(pts[0])
+                     for pts in polygon_points]
+            polygons = list(shapely.polygons(rings))
+        else:
+            polygons = [None] * len(polygon_points)
+            if no_holes:
+                rings = [shapely.linearrings(polygon_points[i][0])
+                         for i in no_holes]
+                batch = shapely.polygons(rings)
+                for idx, poly in zip(no_holes, batch):
+                    polygons[idx] = poly
+            for i, pts in enumerate(polygon_points):
+                if polygons[i] is None:
+                    polygons[i] = Polygon(pts[0], pts[1:])
+    else:
+        # Shapely < 2.0 fallback.
+        polygons = [Polygon(pts[0], pts[1:]) for pts in polygon_points]
 
     df = gpd.GeoDataFrame({column_name: column, "geometry": polygons})
     return df
@@ -823,25 +852,32 @@ def _trace_rings(edge_set):
     return rings
 
 
+@ngjit
 def _simplify_ring(ring):
     """Remove collinear (redundant) vertices from a closed ring."""
     n = len(ring) - 1  # exclude closing point
     if n <= 3:
         return ring
-    keep = []
+    # Single pass into pre-allocated output.
+    out = np.empty((n + 1, 2), dtype=np.float64)
+    count = 0
     for i in range(n):
-        prev = ring[(i - 1) % n]
-        curr = ring[i]
-        nxt = ring[(i + 1) % n]
-        if not ((prev[0] == curr[0] == nxt[0]) or
-                (prev[1] == curr[1] == nxt[1])):
-            keep.append(curr)
-    if len(keep) < n:
-        keep.append(keep[0])
-        return np.array(keep, dtype=np.float64)
+        pi = (i - 1) % n
+        ni = (i + 1) % n
+        if not ((ring[pi, 0] == ring[i, 0] == ring[ni, 0]) or
+                (ring[pi, 1] == ring[i, 1] == ring[ni, 1])):
+            out[count, 0] = ring[i, 0]
+            out[count, 1] = ring[i, 1]
+            count += 1
+    if count < n:
+        out[count, 0] = out[0, 0]
+        out[count, 1] = out[0, 1]
+        count += 1
+        return out[:count].copy()
     return ring
 
 
+@ngjit
 def _signed_ring_area(ring):
     """Shoelace formula for signed area of a closed ring."""
     x = ring[:, 0]
@@ -849,6 +885,7 @@ def _signed_ring_area(ring):
     return 0.5 * (np.dot(x[:-1], y[1:]) - np.dot(x[1:], y[:-1]))
 
 
+@ngjit
 def _point_in_ring(px, py, ring):
     """Ray-casting point-in-polygon test."""
     n = len(ring) - 1

From c770b6ea6493070927549c1fe463690005c22657 Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Sun, 15 Mar 2026 14:39:13 -0700
Subject: [PATCH 2/3] Add performance regression tests for polygonize (#1008)

- Buffer growth: snake-shaped polygon with >64 boundary points
- JIT merge helpers: direct tests of _simplify_ring, _signed_ring_area,
  _point_in_ring
- Dask merge: checkerboard pattern forcing many boundary merges
- Geopandas batch: mixed hole-free and holed polygons through the
  shapely.polygons() batch path
---
 xrspatial/tests/test_polygonize.py | 112 +++++++++++++++++++++++++++++
 1 file changed, 112 insertions(+)

diff --git a/xrspatial/tests/test_polygonize.py b/xrspatial/tests/test_polygonize.py
index 0d581257..1feee493 100644
--- a/xrspatial/tests/test_polygonize.py
+++ b/xrspatial/tests/test_polygonize.py
@@ -457,3 +457,115 @@ def test_polygonize_dask_cupy_matches_numpy(chunks):
         assert val in areas_dcp, f"Value {val} missing from dask+cupy result"
         assert_allclose(areas_dcp[val], areas_np[val],
                         err_msg=f"Area mismatch for value {val}")
+
+
+# --- Performance-related regression tests (#1008) ---
+
+def test_polygonize_1008_large_boundary_buffer_growth():
+    """Single-pass _follow buffer growth: polygon with > 64 boundary points.
+
+    A thin snake-like region forces the boundary tracer to produce many
+    vertices, exercising the dynamic buffer doubling in _follow.
+    """
+    # Horizontal snake: 1-pixel-wide path zigzagging across a 60-column raster.
+    data = np.zeros((6, 60), dtype=np.int32)
+    data[0, :] = 1          # row 0: left to right
+    data[1, 59] = 1         # turn down
+    data[2, :] = 1          # row 2: right to left (fills whole row)
+    data[3, 0] = 1          # turn down
+    data[4, :] = 1          # row 4: left to right
+
+    raster = xr.DataArray(data)
+    values, polygons = polygonize(raster, connectivity=4)
+
+    # The value-1 polygon should have many boundary points (> 64).
+    val1_idx = [i for i, v in enumerate(values) if v == 1]
+    assert len(val1_idx) >= 1
+    for idx in val1_idx:
+        for ring in polygons[idx]:
+            assert ring.shape[1] == 2
+            assert np.array_equal(ring[0], ring[-1])
+
+    # Total area must equal raster size.
+    total_area = sum(
+        assert_polygon_valid_and_get_area(p) for p in polygons)
+    assert_allclose(total_area, data.size)
+
+
+def test_polygonize_1008_jit_merge_helpers():
+    """JIT-compiled _simplify_ring, _signed_ring_area, _point_in_ring."""
+    from ..polygonize import _point_in_ring, _signed_ring_area, _simplify_ring
+
+    # Unit square: CCW exterior.
+    square = np.array([
+        [0, 0], [1, 0], [1, 1], [0, 1], [0, 0]], dtype=np.float64)
+
+    assert_allclose(_signed_ring_area(square), 1.0)
+
+    # Point inside.
+    assert _point_in_ring(0.5, 0.5, square) is True
+    # Point outside.
+    assert _point_in_ring(2.0, 2.0, square) is False
+
+    # Square with collinear midpoints on each edge.
+    with_collinear = np.array([
+        [0, 0], [0.5, 0], [1, 0], [1, 0.5], [1, 1],
+        [0.5, 1], [0, 1], [0, 0.5], [0, 0]], dtype=np.float64)
+    simplified = _simplify_ring(with_collinear)
+    # Should remove the midpoints, leaving 4 corners + closing point.
+    assert simplified.shape == (5, 2)
+    assert_allclose(_signed_ring_area(simplified), 1.0)
+
+    # Ring with no collinear points should be returned unchanged.
+    triangle = np.array([
+        [0, 0], [2, 0], [1, 2], [0, 0]], dtype=np.float64)
+    assert _simplify_ring(triangle) is triangle
+
+
+@dask_array_available
+def test_polygonize_1008_dask_merge_many_boundary_polygons():
+    """Dask merge with many boundary-crossing polygons of the same value.
+
+    Checkerboard pattern in small chunks forces many boundary polygons
+    through the merge path, exercising the JIT-compiled helpers.
+    """
+    # 8x8 checkerboard, chunks of 4x4.
+    data = np.zeros((8, 8), dtype=np.int32)
+    data[::2, ::2] = 1
+    data[1::2, 1::2] = 1
+
+    raster_np = xr.DataArray(data)
+    vals_np, polys_np = polygonize(raster_np, connectivity=4)
+    areas_np = _area_by_value(vals_np, polys_np)
+
+    raster_da = xr.DataArray(da.from_array(data, chunks=(4, 4)))
+    vals_da, polys_da = polygonize(raster_da, connectivity=4)
+    areas_da = _area_by_value(vals_da, polys_da)
+
+    for val in areas_np:
+        assert val in areas_da
+        assert_allclose(areas_da[val], areas_np[val],
+                        err_msg=f"Area mismatch for value {val}")
+
+
+@pytest.mark.skipif(gpd is None, reason="geopandas not installed")
+def test_polygonize_1008_geopandas_batch_with_holes():
+    """Batch shapely construction: mix of hole-free and holed polygons."""
+    # Outer ring of 0s with inner block of 1s containing a 2 (hole in 1).
+    data = np.zeros((6, 6), dtype=np.int32)
+    data[1:5, 1:5] = 1
+    data[2:4, 2:4] = 2
+
+    raster = xr.DataArray(data)
+    df = polygonize(raster, return_type="geopandas", connectivity=4)
+
+    assert isinstance(df, gpd.GeoDataFrame)
+    assert len(df) == 3  # values 0, 1, 2
+
+    # Value 1 polygon should have a hole (the 2-region).
+    row_1 = df[df.DN == 1].iloc[0]
+    geom = row_1.geometry
+    assert len(list(geom.interiors)) == 1
+
+    # Total area should equal raster size.
+    assert_allclose(df.geometry.area.sum(), 36.0)

From 9e4f01cfff7187cf671c7a7b6a01be7f8e268908 Mon Sep 17 00:00:00 2001
From: Brendan Collins <brendancol@gmail.com>
Date: Sun, 15 Mar 2026 18:23:07 -0700
Subject: [PATCH 3/3] Revert single-pass _follow, keep JIT merge helpers and
 batch shapely (#1008)

Benchmarks showed the single-pass _follow was 15-30% slower than the
original two-pass version. The buffer growth check in the inner loop
adds overhead that numba doesn't optimize away, and the two-pass
approach benefits from the data being in cache on the second pass.

Reverted _follow to the original two-pass implementation. The JIT
merge helpers (2.3-2.6x dask speedup) and batch shapely construction
(1.3-1.6x geopandas speedup) are kept.
---
 xrspatial/polygonize.py            | 188 ++++++++++++++---------------
 xrspatial/tests/test_polygonize.py |  31 -----
 2 files changed, 91 insertions(+), 128 deletions(-)

diff --git a/xrspatial/polygonize.py b/xrspatial/polygonize.py
index 439db993..3286d542 100644
--- a/xrspatial/polygonize.py
+++ b/xrspatial/polygonize.py
@@ -107,10 +107,10 @@ def _min_and_max(value0, value1):
 # facing W, and if hole is False the start is on the S edge of the pixel
 # facing E.
 #
-# Single pass with a dynamically-grown buffer (doubling strategy).  Buffer
-# starts at 64 points and grows as needed.  Amortized cost is O(1) per point
-# from O(log n) reallocations, which is cheaper than retracing the boundary
-# a second time.
+# There are two passes.  First pass determines the number of points and
+# allocates a numpy array big enough to take the points, the second pass saves
+# the points.  This is better than a single pass that repeatedly reallocates
+# space for the increasing number of points.
 #
 # Returns the region ID and the 2D array of boundary points.  The last
 # boundary point is the same as the first.
@@ -125,104 +125,98 @@ def _follow(
 ) -> Tuple[int, np.ndarray]:
     region = regions[ij]
     n = nx*ny
+    points = None  # Declare before loop for numba
 
-    # Dynamic buffer for point coordinates (x, y pairs stored flat).
-    buf_cap = 64
-    points = np.empty(2 * buf_cap)
-    npoints = 0
-
-    prev_forward = 0  # Start with an invalid direction.
-    if hole:
-        forward = -1  # Facing W along N edge.
-        left = -nx
-    else:
-        forward = 1  # Facing E along S edge.
-        left = nx
-
-    start_forward = forward
-    start_ij = ij
-
-    while True:
-        if forward == 1 and not hole:
-            # Mark pixel so that it will not be considered a future
-            # non-hole starting pixel.
-            visited[ij] |= 1
-        elif forward == -1 and ij+nx < n:
-            # Mark pixel so that it will not be considered a future
-            # hole starting pixel.
-            visited[ij+nx] |= 2
-
-        if prev_forward != forward:
-            # Grow buffer if needed.
-            if npoints >= buf_cap:
-                old_points = points
-                buf_cap *= 2
-                points = np.empty(2 * buf_cap)
-                points[:2*npoints] = old_points[:2*npoints]
-
-            # Add point.
-            i = ij % nx
-            j = ij // nx
-            if forward == -1:
-                i += 1
-                j += 1
-            elif forward == nx:
-                i += 1
-            elif forward == -nx:
-                j += 1
-            points[2*npoints] = i
-            points[2*npoints+1] = j
-            npoints += 1
-
-        prev_forward = forward
-        ijnext = ij + forward
-        ijnext_right = ijnext - left
-
-        # Determine direction of turn.
-        if abs(forward) == 1:  # Facing E (forward == 1) or W (forward -1)
-            if _diff_row(ij, ijnext, nx):
-                turn = Turn.Left
-            elif (not _outside_domain(ijnext_right, n) and
-                    regions[ijnext_right] == region):
-                turn = Turn.Right
-            elif regions[ijnext] == region:
-                turn = Turn.Straight
-            else:
-                turn = Turn.Left
-        else:  # Facing N (forward == nx) or S (forward == -nx)
-            if _outside_domain(ijnext, n):
-                turn = Turn.Left
-            elif (not _diff_row(ijnext, ijnext_right, nx) and
-                    regions[ijnext_right] == region):
-                turn = Turn.Right
-            elif regions[ijnext] == region:
-                turn = Turn.Straight
-            else:
-                turn = Turn.Left
+    for pass_ in range(2):
+        prev_forward = 0  # Start with an invalid direction.
+        if hole:
+            forward = -1  # Facing W along N edge.
+            left = -nx
+        else:
+            forward = 1  # Facing E along S edge.
+            left = nx
+
+        start_forward = forward
+        start_ij = ij
+        npoints = 0
+
+        while True:
+            if pass_ == 1:
+                if forward == 1 and not hole:
+                    # Mark pixel so that it will not be considered a future
+                    # non-hole starting pixel.
+                    visited[ij] |= 1
+                elif forward == -1 and ij+nx < n:
+                    # Mark pixel so that is will not be considered a future
+                    # hole starting pixel.
+                    visited[ij+nx] |= 2
+
+            if prev_forward != forward:
+                if pass_ == 1:
+                    # Add point.
+                    i = ij % nx
+                    j = ij // nx
+                    if forward == -1:
+                        i += 1
+                        j += 1
+                    elif forward == nx:
+                        i += 1
+                    elif forward == -nx:
+                        j += 1
+                    points[2*npoints] = i
+                    points[2*npoints+1] = j
+                npoints += 1
 
-        # Apply turn.
-        if turn == Turn.Straight:
-            ij = ijnext
-        elif turn == Turn.Left:
-            prev_forward = forward
-            forward = left
-            left = -prev_forward
-        else:  # Turn.Right
             prev_forward = forward
-            forward = -left
-            left = prev_forward
-            ij = ijnext_right
+            ijnext = ij + forward
+            ijnext_right = ijnext - left
+
+            # Determine direction of turn.
+            if abs(forward) == 1:  # Facing E (forward == 1) or W (forward -1)
+                if _diff_row(ij, ijnext, nx):
+                    turn = Turn.Left
+                elif (not _outside_domain(ijnext_right, n) and
+                        regions[ijnext_right] == region):
+                    turn = Turn.Right
+                elif regions[ijnext] == region:
+                    turn = Turn.Straight
+                else:
+                    turn = Turn.Left
+            else:  # Facing N (forward == nx) or S (forward == -nx)
+                if _outside_domain(ijnext, n):
+                    turn = Turn.Left
+                elif (not _diff_row(ijnext, ijnext_right, nx) and
+                        regions[ijnext_right] == region):
+                    turn = Turn.Right
+                elif regions[ijnext] == region:
+                    turn = Turn.Straight
+                else:
+                    turn = Turn.Left
+
+            # Apply turn.
+            if turn == Turn.Straight:
+                ij = ijnext
+            elif turn == Turn.Left:
+                prev_forward = forward
+                forward = left
+                left = -prev_forward
+            else:  # Turn.Right
+                prev_forward = forward
+                forward = -left
+                left = prev_forward
+                ij = ijnext_right
+
+            # Finished boundary when returned to start.
+            if ij == start_ij and forward == start_forward:
+                break
 
-        # Finished boundary when returned to start.
-        if ij == start_ij and forward == start_forward:
-            break
+        if pass_ == 0:
+            # End of first pass.
+            points = np.empty(2*(npoints+1))  # Note extra point at end.
 
-    # Trim buffer to actual size and add closing point.
-    result = np.empty(2*(npoints+1))
-    result[:2*npoints] = points[:2*npoints]
-    result = result.reshape((-1, 2))
-    result[-1] = result[0]  # End point the same as start point.
-    return region, result
+    points = points.reshape((-1, 2))
+    points[-1] = points[0]  # End point the same as start point.
+    return region, points
 
 
 # Generator of numba-compatible comparison functions for values.
diff --git a/xrspatial/tests/test_polygonize.py b/xrspatial/tests/test_polygonize.py
index 1feee493..bf27a784 100644
--- a/xrspatial/tests/test_polygonize.py
+++ b/xrspatial/tests/test_polygonize.py
@@ -461,37 +461,6 @@ def test_polygonize_dask_cupy_matches_numpy(chunks):
 
 # --- Performance-related regression tests (#1008) ---
 
-def test_polygonize_1008_large_boundary_buffer_growth():
-    """Single-pass _follow buffer growth: polygon with > 64 boundary points.
-
-    A thin snake-like region forces the boundary tracer to produce many
-    vertices, exercising the dynamic buffer doubling in _follow.
-    """
-    # Horizontal snake: 1-pixel-wide path zigzagging across a 60-column raster.
-    data = np.zeros((6, 60), dtype=np.int32)
-    data[0, :] = 1          # row 0: left to right
-    data[1, 59] = 1         # turn down
-    data[2, :] = 1          # row 2: right to left (fills whole row)
-    data[3, 0] = 1          # turn down
-    data[4, :] = 1          # row 4: left to right
-
-    raster = xr.DataArray(data)
-    values, polygons = polygonize(raster, connectivity=4)
-
-    # The value-1 polygon should have many boundary points (> 64).
-    val1_idx = [i for i, v in enumerate(values) if v == 1]
-    assert len(val1_idx) >= 1
-    for idx in val1_idx:
-        for ring in polygons[idx]:
-            assert ring.shape[1] == 2
-            assert np.array_equal(ring[0], ring[-1])
-
-    # Total area must equal raster size.
-    total_area = sum(
-        assert_polygon_valid_and_get_area(p) for p in polygons)
-    assert_allclose(total_area, data.size)
-
-
 def test_polygonize_1008_jit_merge_helpers():
     """JIT-compiled _simplify_ring, _signed_ring_area, _point_in_ring."""
     from ..polygonize import _point_in_ring, _signed_ring_area, _simplify_ring