rasterize/kde/line_density: set output attrs['res'] when grid resolution changes (#3571)

xrspatial/kde.py

@@ -546,7 +546,9 @@ def kde(
     Returns
     -------
     xr.DataArray
-        2-D density surface.
+        2-D density surface.  When *template* is omitted, ``attrs['res']``
+        carries the ``(x, y)`` cell size of the grid built from the data
+        extent; when *template* is given, the output inherits its attrs.
     """
     # -- Resolve GeoDataFrame input ----------------------------------------
     if is_geodataframe(x):
@@ -669,10 +671,15 @@ def kde(
             coords=template.coords, dims=template.dims,
             attrs=template.attrs,
         )
+    # No template: the grid was built from x_range/y_range above, so the
+    # cell spacing `dx`/`dy` is the output resolution.  Record it as
+    # `res` so downstream tools (which prefer `attrs['res']` over deriving
+    # cellsize from coords) read the true spacing.
     return xr.DataArray(
         data, name=name,
         dims=['y', 'x'],
         coords={'y': y_coords, 'x': x_coords},
+        attrs={'res': (abs(dx), abs(dy))},
     )
 
 
@@ -725,7 +732,10 @@ def line_density(
     Returns
     -------
     xr.DataArray
-        2-D line-density surface.
+        2-D line-density surface.  When *template* is omitted,
+        ``attrs['res']`` carries the ``(x, y)`` cell size of the grid built
+        from the data extent; when *template* is given, the output inherits
+        its attrs.
     """
     x1a = np.asarray(x1, dtype=np.float64).ravel()
     y1a = np.asarray(y1, dtype=np.float64).ravel()
@@ -799,10 +809,13 @@ def line_density(
             coords=template.coords, dims=template.dims,
             attrs=template.attrs,
         )
+    # No template: `dx`/`dy` is the output cell spacing, so record it as
+    # `res` for downstream tools that prefer `attrs['res']` over coords.
     return xr.DataArray(
         out, name=name,
         dims=['y', 'x'],
         coords={'y': y_coords, 'x': x_coords},
+        attrs={'res': (abs(dx), abs(dy))},
     )
 
 

xrspatial/rasterize.py

@@ -4354,7 +4354,11 @@ def rasterize(
         GeoDataFrame with a ``.crs`` and neither ``like`` nor its
         ``spatial_ref`` coord supplies a CRS, the geometry CRS is
         emitted on the output (``attrs['crs']`` as an EPSG int when one
-        resolves, else ``attrs['crs_wkt']`` as WKT).
+        resolves, else ``attrs['crs_wkt']`` as WKT).  ``attrs['res']``
+        carries the ``(x, y)`` cell size of the output grid.  When the
+        output reuses ``like``'s grid bit-identically the template's
+        ``res`` and ``transform`` are kept; otherwise ``res`` reflects the
+        freshly-built grid and the stale ``transform`` is dropped.
 
     Examples
     --------
@@ -4837,18 +4841,21 @@ def rasterize(
     out_attrs = like_attrs if like_attrs is not None else {}
     for k in ('nodata', '_FillValue', 'nodatavals'):
         out_attrs.pop(k, None)
-    # Grid-shape attrs (res, transform) describe the template's grid.
-    # When the caller overrides bounds/width/height/resolution so the
-    # output grid no longer matches ``like``, leaving the inherited
-    # values in place lies to downstream consumers: get_dataarray_
-    # resolution() prefers ``attrs['res']`` over computing from coords,
-    # so a stale res silently poisons slope/aspect/proximity callers
-    # with the template's cellsize instead of the actual one.  Drop them
-    # when the grid was reshaped; keep them when the output reuses the
-    # template's coords bit-identically.
-    if like_attrs is not None and not reuse_like_coords:
-        for k in ('res', 'transform'):
-            out_attrs.pop(k, None)
+    # Grid-shape attrs (res, transform) describe the output grid.  When
+    # the output reuses ``like``'s coords bit-identically, the template's
+    # res and transform still apply, so leave them alone.  Otherwise the
+    # grid was just built fresh above (the caller overrode
+    # bounds/width/height/resolution, or rasterized without ``like`` at
+    # all), so write the actual cellsize: get_dataarray_resolution()
+    # prefers ``attrs['res']`` over computing from coords, so a stale
+    # template res -- or no res at all -- would feed slope/aspect/
+    # proximity callers the wrong cellsize.  ``px``/``py`` come from the
+    # reshape branch above and are defined whenever ``not
+    # reuse_like_coords``.  Drop any inherited ``transform`` since its
+    # origin/scale no longer match the new grid.
+    if not reuse_like_coords:
+        out_attrs['res'] = (abs(px), abs(py))
+        out_attrs.pop('transform', None)
     try:
         fill_as_float = float(fill)
         fill_is_nan = np.isnan(fill_as_float)

xrspatial/tests/test_kde.py

@@ -497,3 +497,68 @@ def test_cupy_matches_numpy(self, point_cluster, simple_grid, kernel):
         # Compact kernels match exactly.
         tol = 1e-2 if kernel == 'gaussian' else 1e-6
         np.testing.assert_allclose(result_np, np_result.values, rtol=tol)
+
+
+# ---------------------------------------------------------------------------
+# Output resolution metadata (issue #3571)
+# ---------------------------------------------------------------------------
+
+class TestKDEResolutionAttr:
+    """When ``template`` is omitted, kde/line_density build a fresh grid
+    from ``x_range``/``y_range`` and must record its cell spacing as
+    ``attrs['res']`` so downstream tools (which prefer ``attrs['res']``
+    over deriving cellsize from coords) read the true resolution.
+    """
+
+    def test_kde_sets_res_from_range(self):
+        # x spans 0..100 over 11 cells -> dx = 10; y spans 0..50 over 6
+        # cells -> dy = 10.  Rectangular grid pins the (x_res, y_res) order.
+        x = np.array([10.0, 20.0, 80.0])
+        y = np.array([5.0, 10.0, 40.0])
+        result = kde(x, y, bandwidth=5.0,
+                     x_range=(0, 100), y_range=(0, 50),
+                     width=11, height=6)
+        assert result.attrs['res'] == (10.0, 10.0)
+
+    def test_kde_res_matches_coords(self):
+        from xrspatial.utils import get_dataarray_resolution
+        x = np.array([1.0, 2.0, 3.0])
+        y = np.array([1.0, 2.0, 3.0])
+        result = kde(x, y, bandwidth=1.0,
+                     x_range=(0, 20), y_range=(0, 10),
+                     width=21, height=6)
+        # res must agree with the cellsize implied by the coordinates.
+        cx = float(result.x.values[1] - result.x.values[0])
+        cy = float(result.y.values[1] - result.y.values[0])
+        assert result.attrs['res'] == (abs(cx), abs(cy))
+        # get_dataarray_resolution reads res directly.
+        rx, ry = get_dataarray_resolution(result)
+        assert (rx, ry) == result.attrs['res']
+
+    def test_kde_template_res_from_template(self, simple_grid):
+        # A template carrying res propagates it (output grid == template).
+        template = simple_grid.copy()
+        template.attrs['res'] = (0.5, 0.5)
+        x = np.array([0.0, 1.0])
+        y = np.array([0.0, 1.0])
+        result = kde(x, y, bandwidth=1.0, template=template)
+        assert result.attrs['res'] == (0.5, 0.5)
+
+    def test_line_density_sets_res_from_range(self):
+        result = line_density([0.0], [0.0], [100.0], [50.0],
+                              bandwidth=5.0,
+                              x_range=(0, 100), y_range=(0, 50),
+                              width=11, height=6)
+        assert result.attrs['res'] == (10.0, 10.0)
+
+    def test_line_density_res_matches_coords(self):
+        from xrspatial.utils import get_dataarray_resolution
+        result = line_density([0.0], [0.0], [20.0], [10.0],
+                              bandwidth=2.0,
+                              x_range=(0, 20), y_range=(0, 10),
+                              width=21, height=6)
+        cx = float(result.x.values[1] - result.x.values[0])
+        cy = float(result.y.values[1] - result.y.values[0])
+        assert result.attrs['res'] == (abs(cx), abs(cy))
+        rx, ry = get_dataarray_resolution(result)
+        assert (rx, ry) == result.attrs['res']

xrspatial/tests/test_rasterize.py

@@ -1913,12 +1913,15 @@ def test_fill_value_omitted_for_nan(self):
         assert 'nodatavals' not in result.attrs
 
     def test_no_like_no_attrs_pollution(self):
-        """Without `like` and with NaN fill, attrs must stay empty."""
+        """Without `like` and with NaN fill, the only attr is the grid
+        ``res`` describing the freshly-built output; no nodata keys leak.
+        """
         result = rasterize(
             [(box(2, 2, 8, 8), 1.0)],
             width=10, height=10, bounds=(0, 0, 10, 10),
         )
-        assert result.attrs == {}
+        # 10x10 over a 10-unit span -> 1.0/pixel.
+        assert result.attrs == {'res': (1.0, 1.0)}
 
     @skip_no_dask
     def test_like_attrs_propagated_dask(self):
@@ -2047,12 +2050,13 @@ def test_geotiff_round_trip_preserves_fill(self, tmp_path):
 
 
 class TestLikeStaleGridAttrs2251:
-    """Issue #2251 -- ``like.attrs['res']`` and ``attrs['transform']`` describe
-    the template's grid.  When the caller overrides ``bounds``,
+    """Issue #2251 / #3571 -- ``like.attrs['res']`` and ``attrs['transform']``
+    describe the template's grid.  When the caller overrides ``bounds``,
     ``width``/``height``, or ``resolution`` so the output grid is no longer
-    identical to ``like``, those keys have to be dropped so downstream
-    consumers (``get_dataarray_resolution`` prefers ``attrs['res']`` over
-    coords) don't see a stale cellsize.
+    identical to ``like``, the template's cellsize no longer applies.
+    ``res`` is refreshed to the grid actually built (so downstream consumers,
+    which prefer ``attrs['res']`` over coords, read the true cellsize) and the
+    inherited ``transform`` is dropped (its origin/scale no longer match).
     """
 
     @staticmethod
@@ -2077,8 +2081,9 @@ def test_bounds_override_strips_stale_grid_attrs(self):
             like=like, bounds=(0, 0, 100, 100),
             width=10, height=10, fill=0,
         )
-        # Grid shape moved from 1.0/pixel to 10.0/pixel; stale attrs gone.
-        assert 'res' not in result.attrs
+        # Grid shape moved from 1.0/pixel to 10.0/pixel; res refreshed to
+        # the new cellsize and the stale transform dropped.
+        assert result.attrs.get('res') == (10.0, 10.0)
         assert 'transform' not in result.attrs
         # Non-grid-shape attrs (crs) still propagate.
         assert result.attrs.get('crs') == 'EPSG:32610'
@@ -2089,7 +2094,8 @@ def test_width_height_override_strips_stale_grid_attrs(self):
             [(box(2, 2, 8, 8), 1.0)],
             like=like, width=5, height=5, fill=0,
         )
-        assert 'res' not in result.attrs
+        # Template spans 0..10 over 5 cells -> 2.0/pixel.
+        assert result.attrs.get('res') == (2.0, 2.0)
         assert 'transform' not in result.attrs
         assert result.attrs.get('crs') == 'EPSG:32610'
 
@@ -2099,7 +2105,7 @@ def test_resolution_override_strips_stale_grid_attrs(self):
             [(box(2, 2, 8, 8), 1.0)],
             like=like, resolution=0.5, fill=0,
         )
-        assert 'res' not in result.attrs
+        assert result.attrs.get('res') == (0.5, 0.5)
         assert 'transform' not in result.attrs
         assert result.attrs.get('crs') == 'EPSG:32610'
 
@@ -2138,8 +2144,8 @@ def test_get_dataarray_resolution_consistent_after_override(self):
             like=like, bounds=(0, 0, 100, 100),
             width=10, height=10, fill=0,
         )
-        # With res stripped, get_dataarray_resolution falls back to
-        # computing from coords; both axes should be 10.0.
+        # res is refreshed to the actual grid, so get_dataarray_resolution
+        # reads it directly; both axes should be 10.0.
         cx, cy = get_dataarray_resolution(result)
         assert abs(cx - 10.0) < 1e-9
         assert abs(cy - 10.0) < 1e-9
@@ -2152,7 +2158,7 @@ def test_bounds_override_strips_stale_grid_attrs_dask(self):
             like=like, bounds=(0, 0, 100, 100),
             width=10, height=10, fill=0, chunks=5,
         )
-        assert 'res' not in result.attrs
+        assert result.attrs.get('res') == (10.0, 10.0)
         assert 'transform' not in result.attrs
         assert result.attrs.get('crs') == 'EPSG:32610'
 
@@ -2164,7 +2170,7 @@ def test_bounds_override_strips_stale_grid_attrs_cupy(self):
             like=like, bounds=(0, 0, 100, 100),
             width=10, height=10, fill=0, gpu=True,
         )
-        assert 'res' not in result.attrs
+        assert result.attrs.get('res') == (10.0, 10.0)
         assert 'transform' not in result.attrs
         assert result.attrs.get('crs') == 'EPSG:32610'
 
@@ -2177,7 +2183,7 @@ def test_bounds_override_strips_stale_grid_attrs_dask_cupy(self):
             like=like, bounds=(0, 0, 100, 100),
             width=10, height=10, fill=0, gpu=True, chunks=5,
         )
-        assert 'res' not in result.attrs
+        assert result.attrs.get('res') == (10.0, 10.0)
         assert 'transform' not in result.attrs
         assert result.attrs.get('crs') == 'EPSG:32610'
 

xrspatial/tests/test_rasterize_coverage_2026_05_21.py

@@ -483,12 +483,11 @@ def test_rectangular_pixels_polygon_parity(self, backend_name, kw):
             err_msg=f"{backend_name} disagreed on (rx=2, ry=1) cell size")
 
     def test_rectangular_pixels_attrs_res(self):
-        """The resolved cell size lands on the output if a like attrs
-        dict is supplied with res; without like, no res attr is set --
-        this test pins that contract so callers know which path emits
-        the metadata."""
-        # No like -> no res attr.
+        """The resolved cell size always lands on the output ``res`` attr,
+        even without ``like`` (issue #3571).  Rectangular pixels also pin
+        the (x_res, y_res) ordering: resolution=(2, 1) over (0,0)-(10,10)
+        builds a 5x10 grid, so res is (2.0, 1.0)."""
         r = rasterize([(box(2, 2, 8, 8), 3.0)],
                       resolution=(2.0, 1.0),
                       bounds=(0, 0, 10, 10), fill=0)
-        assert 'res' not in r.attrs
+        assert r.attrs['res'] == (2.0, 1.0)