Recompute res attribute for downsampled preview() output

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

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))