Fix KDE all-zero output with descending-coordinate templates (#1199)

* 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 KDE all-zero output with descending-coordinate templates (#1198)

The bounding-box index calculation in _kde_cpu and _line_density_cpu
divided by dx/dy to convert coordinate offsets to pixel indices.
When dy or dx was negative (descending coordinates, common for
north-up rasters), the division flipped lo/hi so the inner loops
never executed, producing all-zero output.

Fix: compute both index endpoints and use min/max to get the correct
lo and hi regardless of spacing sign.

Author: brendancol

xrspatial/kde.py

@@ -118,10 +118,16 @@ def _kde_cpu(xs, ys, ws, out, x0, y0, dx, dy, bw, kernel_id):
         w = ws[p]
 
         # Pixel range that falls within the cutoff.
-        col_lo = int(max(0, (px - cutoff - x0) / dx))
-        col_hi = int(min(cols - 1, (px + cutoff - x0) / dx)) + 1
-        row_lo = int(max(0, (py - cutoff - y0) / dy))
-        row_hi = int(min(rows - 1, (py + cutoff - y0) / dy)) + 1
+        # Compute both endpoints and use min/max so that negative
+        # spacing (descending coordinates) still produces lo <= hi.
+        c_a = int((px - cutoff - x0) / dx)
+        c_b = int((px + cutoff - x0) / dx)
+        col_lo = max(0, min(c_a, c_b))
+        col_hi = min(cols - 1, max(c_a, c_b)) + 1
+        r_a = int((py - cutoff - y0) / dy)
+        r_b = int((py + cutoff - y0) / dy)
+        row_lo = max(0, min(r_a, r_b))
+        row_hi = min(rows - 1, max(r_a, r_b)) + 1
 
         for r in range(row_lo, row_hi):
             cy = y0 + r * dy
@@ -200,10 +206,14 @@ def _line_density_cpu(x1s, y1s, x2s, y2s, ws, out,
             px = ax + t * (bx - ax)
             py = ay + t * (by - ay)
 
-            col_lo = int(max(0, (px - cutoff - x0) / dx))
-            col_hi = int(min(cols - 1, (px + cutoff - x0) / dx)) + 1
-            row_lo = int(max(0, (py - cutoff - y0) / dy))
-            row_hi = int(min(rows - 1, (py + cutoff - y0) / dy)) + 1
+            c_a = int((px - cutoff - x0) / dx)
+            c_b = int((px + cutoff - x0) / dx)
+            col_lo = max(0, min(c_a, c_b))
+            col_hi = min(cols - 1, max(c_a, c_b)) + 1
+            r_a = int((py - cutoff - y0) / dy)
+            r_b = int((py + cutoff - y0) / dy)
+            row_lo = max(0, min(r_a, r_b))
+            row_hi = min(rows - 1, max(r_a, r_b)) + 1
 
             for r in range(row_lo, row_hi):
                 cy = y0 + r * dy

xrspatial/tests/test_kde.py

@@ -198,6 +198,163 @@ def test_template_not_2d_raises(self):
             kde([0], [0], bandwidth=1.0, template=t)
 
 
+# ---------------------------------------------------------------------------
+# Descending coordinate templates (#1198)
+# ---------------------------------------------------------------------------
+
+class TestDescendingCoordinates:
+    """KDE must produce correct results when template coords are descending.
+
+    Geospatial rasters commonly have row 0 at the top (descending y).
+    Before the fix in #1198, negative dy/dx caused the bounding-box
+    index math to produce lo > hi, so the inner loops never ran.
+    """
+
+    @pytest.fixture
+    def asc_template(self):
+        return xr.DataArray(
+            np.zeros((16, 16), dtype=np.float64),
+            dims=['y', 'x'],
+            coords={
+                'y': np.linspace(-4, 4, 16),
+                'x': np.linspace(-4, 4, 16),
+            },
+        )
+
+    @pytest.fixture
+    def desc_y_template(self):
+        """Row 0 = top (descending y, ascending x)."""
+        return xr.DataArray(
+            np.zeros((16, 16), dtype=np.float64),
+            dims=['y', 'x'],
+            coords={
+                'y': np.linspace(4, -4, 16),
+                'x': np.linspace(-4, 4, 16),
+            },
+        )
+
+    @pytest.fixture
+    def desc_x_template(self):
+        """Ascending y, descending x."""
+        return xr.DataArray(
+            np.zeros((16, 16), dtype=np.float64),
+            dims=['y', 'x'],
+            coords={
+                'y': np.linspace(-4, 4, 16),
+                'x': np.linspace(4, -4, 16),
+            },
+        )
+
+    @pytest.fixture
+    def desc_both_template(self):
+        """Both axes descending."""
+        return xr.DataArray(
+            np.zeros((16, 16), dtype=np.float64),
+            dims=['y', 'x'],
+            coords={
+                'y': np.linspace(4, -4, 16),
+                'x': np.linspace(4, -4, 16),
+            },
+        )
+
+    def test_descending_y_nonzero(self, desc_y_template):
+        result = kde([0.0], [0.0], bandwidth=1.0, template=desc_y_template)
+        assert float(result.sum()) > 0.0
+
+    def test_descending_y_matches_ascending(self, asc_template, desc_y_template):
+        """Descending-y result should equal ascending-y result, row-flipped."""
+        r_asc = kde([0.0], [0.0], bandwidth=1.0, template=asc_template)
+        r_desc = kde([0.0], [0.0], bandwidth=1.0, template=desc_y_template)
+        np.testing.assert_allclose(
+            r_desc.values[::-1], r_asc.values, rtol=1e-12,
+        )
+
+    def test_descending_x_matches_ascending(self, asc_template, desc_x_template):
+        """Descending-x result should equal ascending-x result, col-flipped."""
+        r_asc = kde([0.0], [0.0], bandwidth=1.0, template=asc_template)
+        r_desc = kde([0.0], [0.0], bandwidth=1.0, template=desc_x_template)
+        np.testing.assert_allclose(
+            r_desc.values[:, ::-1], r_asc.values, rtol=1e-12,
+        )
+
+    def test_both_descending_matches_ascending(self, asc_template,
+                                                desc_both_template):
+        r_asc = kde([0.0], [0.0], bandwidth=1.0, template=asc_template)
+        r_both = kde([0.0], [0.0], bandwidth=1.0, template=desc_both_template)
+        np.testing.assert_allclose(
+            r_both.values[::-1, ::-1], r_asc.values, rtol=1e-12,
+        )
+
+    @pytest.mark.parametrize('kernel', ['epanechnikov', 'quartic'])
+    def test_compact_kernels_descending_y(self, asc_template,
+                                           desc_y_template, kernel):
+        """Compact kernels match exactly (no cutoff-boundary rounding)."""
+        pts_x = [0.0, 1.0, -1.0]
+        pts_y = [0.0, 1.0, -1.0]
+        r_asc = kde(pts_x, pts_y, bandwidth=1.0, kernel=kernel,
+                    template=asc_template)
+        r_desc = kde(pts_x, pts_y, bandwidth=1.0, kernel=kernel,
+                     template=desc_y_template)
+        np.testing.assert_allclose(
+            r_desc.values[::-1], r_asc.values, rtol=1e-12,
+        )
+
+    def test_gaussian_descending_y_multipoint(self, asc_template,
+                                               desc_y_template):
+        """Gaussian with multiple points: allow tiny diffs at cutoff edge.
+
+        The 4*bw box cutoff means int() truncation can include/exclude
+        a different fringe pixel when spacing flips sign. The affected
+        values are near exp(-8) ~ 3e-4, so atol=1e-5 is plenty tight.
+        """
+        pts_x = [0.0, 1.0, -1.0]
+        pts_y = [0.0, 1.0, -1.0]
+        r_asc = kde(pts_x, pts_y, bandwidth=1.0, kernel='gaussian',
+                    template=asc_template)
+        r_desc = kde(pts_x, pts_y, bandwidth=1.0, kernel='gaussian',
+                     template=desc_y_template)
+        np.testing.assert_allclose(
+            r_desc.values[::-1], r_asc.values, atol=1e-5,
+        )
+
+    def test_line_density_descending_y_quartic(self, asc_template,
+                                                desc_y_template):
+        """Compact kernel: exact match with descending y."""
+        r_asc = line_density([0.0], [0.0], [1.0], [1.0],
+                             bandwidth=0.5, kernel='quartic',
+                             template=asc_template)
+        r_desc = line_density([0.0], [0.0], [1.0], [1.0],
+                              bandwidth=0.5, kernel='quartic',
+                              template=desc_y_template)
+        np.testing.assert_allclose(
+            r_desc.values[::-1], r_asc.values, rtol=1e-12,
+        )
+
+    def test_line_density_descending_y_gaussian(self, asc_template,
+                                                 desc_y_template):
+        """Gaussian: allow tiny diffs at the 4*bw cutoff edge."""
+        r_asc = line_density([0.0], [0.0], [1.0], [1.0],
+                             bandwidth=0.5, kernel='gaussian',
+                             template=asc_template)
+        r_desc = line_density([0.0], [0.0], [1.0], [1.0],
+                              bandwidth=0.5, kernel='gaussian',
+                              template=desc_y_template)
+        np.testing.assert_allclose(
+            r_desc.values[::-1], r_asc.values, atol=1e-4,
+        )
+
+    def test_line_density_descending_x(self, asc_template, desc_x_template):
+        r_asc = line_density([0.0], [0.0], [1.0], [1.0],
+                             bandwidth=0.5, kernel='quartic',
+                             template=asc_template)
+        r_desc = line_density([0.0], [0.0], [1.0], [1.0],
+                              bandwidth=0.5, kernel='quartic',
+                              template=desc_x_template)
+        np.testing.assert_allclose(
+            r_desc.values[:, ::-1], r_asc.values, rtol=1e-12,
+        )
+
+
 # ---------------------------------------------------------------------------
 # Line density
 # ---------------------------------------------------------------------------