Fix crop=True dropping boundary pixels when all_touched=True (#1197) (#1200)

* Fix unbounded allocation DoS and VRT path traversal in geotiff

Two security fixes for the geotiff subpackage:

1. Add a configurable max_pixels guard to read_to_array() and all
   internal read functions (_read_strips, _read_tiles, _read_cog_http).
   A crafted TIFF with fabricated header dimensions could previously
   trigger multi-TB allocations. The default limit is 1 billion pixels
   (~4 GB for float32 single-band), overridable via max_pixels kwarg.
   Fixes #1184.

2. Canonicalize VRT source filenames with os.path.realpath() after
   resolving relative paths. Previously, a VRT file with "../" in
   SourceFilename could read arbitrary files outside the VRT directory.
   Fixes #1185.

* Fix VRT parser test failure on Windows

os.path.realpath() converts Unix-style paths to Windows paths on
Windows (e.g. /data/tile.tif becomes D:\data\tile.tif). Use
os.path.realpath() in the assertion so it matches the production
code's canonicalization on all platforms.

* Fix crop=True dropping boundary pixels when all_touched=True (#1197)

_crop_to_bbox compared pixel center coordinates against the geometry
bounding box without accounting for pixel cell extent. When
all_touched=True, pixels whose centers fell just outside the bbox were
excluded even though their cells overlapped the polygon.

Now _crop_to_bbox receives the all_touched flag and expands the bbox
comparison by half a pixel on each side when set, so rasterize gets to
see every pixel whose cell intersects the geometry.

Also removed dead ascending/descending branches that computed the same
mask regardless of coordinate order.

Author: brendancol

xrspatial/polygon_clip.py

@@ -79,10 +79,14 @@ def _resolve_geometry(geometry):
     )
 
 
-def _crop_to_bbox(raster, geom_pairs):
+def _crop_to_bbox(raster, geom_pairs, all_touched=False):
     """Slice the raster to the bounding box of the geometry.
 
     Returns the sliced DataArray and the geometry pairs (unchanged).
+
+    When *all_touched* is True the bounding-box comparison is expanded by
+    half a pixel on every side so that pixels whose cells overlap the
+    geometry boundary are not prematurely excluded.
     """
     from shapely.ops import unary_union
     merged = unary_union([g for g, _ in geom_pairs])
@@ -91,19 +95,24 @@ def _crop_to_bbox(raster, geom_pairs):
     y_coords = raster.coords[raster.dims[-2]].values
     x_coords = raster.coords[raster.dims[-1]].values
 
-    # Determine coordinate ordering (ascending or descending)
-    y_ascending = y_coords[-1] >= y_coords[0]
-    x_ascending = x_coords[-1] >= x_coords[0]
-
-    if y_ascending:
-        y_mask = (y_coords >= miny) & (y_coords <= maxy)
-    else:
-        y_mask = (y_coords >= miny) & (y_coords <= maxy)
+    # When all_touched is set, expand the bbox by half a pixel so that
+    # pixels whose cells overlap the geometry survive the crop.
+    if all_touched:
+        if len(x_coords) > 1:
+            half_px_x = abs(float(x_coords[1] - x_coords[0])) / 2.0
+        else:
+            half_px_x = 0.0
+        if len(y_coords) > 1:
+            half_py_y = abs(float(y_coords[1] - y_coords[0])) / 2.0
+        else:
+            half_py_y = 0.0
+        minx -= half_px_x
+        maxx += half_px_x
+        miny -= half_py_y
+        maxy += half_py_y
 
-    if x_ascending:
-        x_mask = (x_coords >= minx) & (x_coords <= maxx)
-    else:
-        x_mask = (x_coords >= minx) & (x_coords <= maxx)
+    y_mask = (y_coords >= miny) & (y_coords <= maxy)
+    x_mask = (x_coords >= minx) & (x_coords <= maxx)
 
     y_idx = np.where(y_mask)[0]
     x_idx = np.where(x_mask)[0]
@@ -186,7 +195,7 @@ def clip_polygon(
 
     # Optionally crop to bounding box first (reduces rasterize cost)
     if crop:
-        raster = _crop_to_bbox(raster, geom_pairs)
+        raster = _crop_to_bbox(raster, geom_pairs, all_touched=all_touched)
 
     # Build a binary mask via rasterize, aligned to the (possibly cropped)
     # raster grid.  Always produce a plain numpy mask first, then convert

xrspatial/tests/test_polygon_clip.py

@@ -152,6 +152,42 @@ def test_all_touched(self):
         n_touched = np.count_nonzero(np.isfinite(result_touched.values))
         assert n_touched >= n_default
 
+    def test_all_touched_crop_matches_nocrop(self):
+        """crop=True with all_touched=True must not drop boundary pixels.
+
+        Regression test for #1197: _crop_to_bbox was comparing pixel
+        centers against the geometry bounding box without accounting for
+        pixel cell extent, so pixels whose centers fell just outside the
+        bbox were excluded even though their cells overlapped the polygon.
+        """
+        raster = _make_raster()
+        # Polygon whose edges land between pixel centers on all four
+        # sides.  Pixel spacing is 0.5, so the left edge at x=0.15
+        # sits inside the cell of pixel x=0.0 (cell [-0.25, 0.25]).
+        poly = Polygon([(0.15, 0.15), (0.15, 3.35), (2.35, 3.35),
+                         (2.35, 0.15)])
+
+        result_crop = clip_polygon(raster, poly, crop=True,
+                                   all_touched=True)
+        result_nocrop = clip_polygon(raster, poly, crop=False,
+                                     all_touched=True)
+
+        # The crop path must keep every pixel the nocrop path keeps.
+        n_crop = np.count_nonzero(np.isfinite(result_crop.values))
+        n_nocrop = np.count_nonzero(np.isfinite(result_nocrop.values))
+        assert n_crop == n_nocrop, (
+            f"crop=True lost {n_nocrop - n_crop} boundary pixels"
+        )
+
+        # Pixel values must match wherever both arrays have data.
+        # Align the cropped result back into the full grid for comparison.
+        aligned = result_nocrop.copy()
+        crop_y = result_crop.coords['y'].values
+        crop_x = result_crop.coords['x'].values
+        aligned_slice = aligned.sel(y=crop_y, x=crop_x)
+        np.testing.assert_array_equal(result_crop.values,
+                                      aligned_slice.values)
+
     def test_single_cell_raster(self):
         """1x1 raster with polygon that covers it."""
         data = np.array([[42.0]])
@@ -207,6 +243,23 @@ def test_crop_matches_numpy(self):
             dk_result.values, np_result.values, equal_nan=True
         )
 
+    def test_all_touched_crop_matches_nocrop(self):
+        """Dask: crop=True with all_touched=True keeps boundary pixels (#1197)."""
+        dk_raster = _make_raster(backend='dask+numpy', chunks=(4, 3))
+        poly = Polygon([(0.15, 0.15), (0.15, 3.35), (2.35, 3.35),
+                         (2.35, 0.15)])
+
+        result_crop = clip_polygon(dk_raster, poly, crop=True,
+                                   all_touched=True)
+        result_nocrop = clip_polygon(dk_raster, poly, crop=False,
+                                     all_touched=True)
+
+        n_crop = np.count_nonzero(np.isfinite(result_crop.values))
+        n_nocrop = np.count_nonzero(np.isfinite(result_nocrop.values))
+        assert n_crop == n_nocrop, (
+            f"crop=True lost {n_nocrop - n_crop} boundary pixels"
+        )
+
 
 # ---------------------------------------------------------------------------
 # CuPy backend