From d53cf61f577a9f7a04617d85d5b8268d02de8a0e Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 28 Jun 2026 20:47:56 -0700 Subject: [PATCH 1/2] Update output attrs['res'] when resolution changes in rasterize/kde (#3571) rasterize(): the reshape path (caller overrides bounds/width/height/ resolution, or rasterizes without like=) now records the freshly-built grid's cell size in attrs['res'] instead of leaving it absent. The reuse-like-coords path still keeps the template's res/transform. Stale inherited transform is still dropped on reshape. kde()/line_density(): when template is omitted, set attrs['res'] to the (dx, dy) cell spacing of the grid built from x_range/y_range. get_dataarray_resolution() prefers attrs['res'] over deriving cellsize from coords, so a missing or stale res fed the wrong cellsize to downstream slope/aspect/proximity callers. resample() already did this; this brings rasterize/kde/line_density in line. Updates the issue #2251 tests to assert the refreshed res value rather than its absence. --- xrspatial/kde.py | 17 +++++- xrspatial/rasterize.py | 33 +++++++----- xrspatial/tests/test_kde.py | 53 +++++++++++++++++++ xrspatial/tests/test_rasterize.py | 38 +++++++------ .../test_rasterize_coverage_2026_05_21.py | 11 ++-- 5 files changed, 115 insertions(+), 37 deletions(-) diff --git a/xrspatial/kde.py b/xrspatial/kde.py index ce95ac3ad..ff9e8a2bb 100644 --- a/xrspatial/kde.py +++ b/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))}, ) diff --git a/xrspatial/rasterize.py b/xrspatial/rasterize.py index d001246fb..c6dc3fd43 100644 --- a/xrspatial/rasterize.py +++ b/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) diff --git a/xrspatial/tests/test_kde.py b/xrspatial/tests/test_kde.py index 79c70f1e8..d5041364b 100644 --- a/xrspatial/tests/test_kde.py +++ b/xrspatial/tests/test_kde.py @@ -497,3 +497,56 @@ 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) diff --git a/xrspatial/tests/test_rasterize.py b/xrspatial/tests/test_rasterize.py index baa033b77..dcd4f4254 100644 --- a/xrspatial/tests/test_rasterize.py +++ b/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' diff --git a/xrspatial/tests/test_rasterize_coverage_2026_05_21.py b/xrspatial/tests/test_rasterize_coverage_2026_05_21.py index bcd83c214..5522c70de 100644 --- a/xrspatial/tests/test_rasterize_coverage_2026_05_21.py +++ b/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) From 939d358ae9dba7db5db1fde28c73b190f7f0e905 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 28 Jun 2026 20:51:02 -0700 Subject: [PATCH 2/2] Address review nit: add line_density res-vs-coords test (#3571) --- xrspatial/tests/test_kde.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/xrspatial/tests/test_kde.py b/xrspatial/tests/test_kde.py index d5041364b..9da1c0cd1 100644 --- a/xrspatial/tests/test_kde.py +++ b/xrspatial/tests/test_kde.py @@ -550,3 +550,15 @@ def test_line_density_sets_res_from_range(self): 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']