diff --git a/xrspatial/preview.py b/xrspatial/preview.py index 87b3d9509..74911cac1 100644 --- a/xrspatial/preview.py +++ b/xrspatial/preview.py @@ -282,6 +282,24 @@ def _interpolate_coords(coords, n_out): return np.interp(indices, np.arange(len(vals)), vals.astype(np.float64)) +def _res_from_coords(result, y_dim, x_dim): + """Pixel size ``(res_x, res_y)`` from the output coordinate spacing. + + Returns ``None`` when either axis has fewer than two coordinates, so + the caller can leave any existing ``res`` untouched. Spacing is taken + as a positive magnitude to match how ``res`` is stored elsewhere + (templates, terrain), regardless of whether the axis ascends or + descends. + """ + if x_dim not in result.coords or y_dim not in result.coords: + return None + xs = np.asarray(result.coords[x_dim].values) + ys = np.asarray(result.coords[y_dim].values) + if xs.size < 2 or ys.size < 2: + return None + return (abs(float(xs[1] - xs[0])), abs(float(ys[1] - ys[0]))) + + # --------------------------------------------------------------------------- # Second-pass refinement # --------------------------------------------------------------------------- @@ -350,7 +368,8 @@ def preview(agg, width=1000, height=None, method='mean', name='preview'): Returns ------- xr.DataArray - Downsampled raster with updated coordinates. + Downsampled raster with updated coordinates. ``attrs['res']`` is + recomputed from the output spacing when the input carried a ``res``. """ _validate_raster(agg, func_name='preview', ndim=2) @@ -444,4 +463,12 @@ def preview(agg, width=1000, height=None, method='mean', name='preview'): ) result.name = name + # The output grid is coarser than the input, so the copied 'res' is + # stale. Recompute it from the downsampled coordinates so it matches + # the returned array's x/y spacing. + if 'res' in result.attrs: + new_res = _res_from_coords(result, y_dim, x_dim) + if new_res is not None: + result.attrs = {**result.attrs, 'res': new_res} + return result diff --git a/xrspatial/tests/test_preview.py b/xrspatial/tests/test_preview.py index 876bbe593..4c840bfc1 100644 --- a/xrspatial/tests/test_preview.py +++ b/xrspatial/tests/test_preview.py @@ -480,3 +480,66 @@ def test_accessor_dataset(): result = ds.xrs.preview(width=40) assert isinstance(result, xr.Dataset) assert result['elev'].shape == (20, 40) + + +# ---- res attribute recomputed for the downsampled grid (issue #3569) ---- + +def _expected_res(result): + """Pixel size implied by a result's x/y coordinate spacing.""" + xs = np.asarray(result.coords['x'].values) + ys = np.asarray(result.coords['y'].values) + return (abs(float(xs[1] - xs[0])), abs(float(ys[1] - ys[0]))) + + +@pytest.mark.parametrize( + 'method', ['mean', 'median', 'max', 'min', 'nearest', 'bilinear'], +) +def test_res_matches_downsampled_coords(method): + # create_test_raster seeds res=(0.5, 0.5) with matching 0.5 spacing; + # after a 10x downsample the copied res would be stale. + agg = _make_raster(200, 400) + assert agg.attrs['res'] == (0.5, 0.5) + result = preview(agg, width=40, method=method) + np.testing.assert_allclose(result.attrs['res'], _expected_res(result)) + # the source-grid res must not survive unchanged + assert not np.allclose(result.attrs['res'], (0.5, 0.5)) + + +def test_res_does_not_mutate_input(): + agg = _make_raster(200, 400) + preview(agg, width=40) + assert agg.attrs['res'] == (0.5, 0.5) + + +def test_res_passthrough_unchanged(): + """No downsampling means the original res is returned untouched.""" + agg = _make_raster(5, 5) + result = preview(agg, width=100) + assert result.attrs['res'] == (0.5, 0.5) + + +def test_no_res_attr_stays_absent(): + """preview should not invent a res when the input never had one.""" + data = np.random.default_rng(1).random((200, 400)).astype(np.float32) + agg = create_test_raster(data, attrs={'crs': 'EPSG: 5070'}) + assert 'res' not in agg.attrs + result = preview(agg, width=40) + assert 'res' not in result.attrs + + +@dask_array_available +def test_res_dask_with_refinement(): + # an odd target width forces the dask snap path to overshoot and the + # second-pass refinement to rebuild the coords. + agg = _make_raster(300, 600, backend='dask', chunks=(60, 60)) + result = preview(agg, width=37).compute() + np.testing.assert_allclose(result.attrs['res'], _expected_res(result)) + assert not np.allclose(result.attrs['res'], (0.5, 0.5)) + + +@cuda_and_cupy_available +def test_res_cupy(): + agg = _make_raster(200, 400, backend='cupy') + result = preview(agg, width=40, method='mean') + np.testing.assert_allclose(result.attrs['res'], _expected_res(result)) + assert not np.allclose(result.attrs['res'], (0.5, 0.5))