Honor requested resolution exactly in from_template() (#3494)

xrspatial/templates.py

@@ -182,6 +182,12 @@ def from_template(name, resolution=None, *, preserve=None, backend="numpy",
     (north-up, descending ``y``) and ``res``/``crs`` attributes. It covers the
     study area's rectangular bounding box and is meant as a starting canvas.
 
+    The requested ``resolution`` is honored exactly. The study-area box is
+    rarely a whole number of cells wide, so the far edges (right, top) are
+    nudged out by up to half a cell to land on an exact multiple of the cell
+    size, anchoring the lower-left corner. ``attrs['res']`` therefore matches
+    the ``resolution`` you pass.
+
     Parameters
     ----------
     name : str
@@ -268,10 +274,14 @@ def from_template(name, resolution=None, *, preserve=None, backend="numpy",
             f"Use a coarser resolution."
         )
 
-    ys, xs = _make_output_coords(bounds, (height, width))
-    # Realized cell size from the integer grid (may differ slightly from input).
-    actual_res_x = (right - left) / width
-    actual_res_y = (top - bottom) / height
+    # Honor the requested resolution exactly: anchor the lower-left corner and
+    # nudge the far edges to an exact multiple of the cell size, so res comes
+    # back as (res_x, res_y) instead of drifting when the bbox extent isn't a
+    # whole number of cells. Mirrors the reproject grid path
+    # (_compute_output_grid in reproject/_grid.py).
+    right = left + width * res_x
+    top = bottom + height * res_y
+    ys, xs = _make_output_coords((left, bottom, right, top), (height, width))
 
     data = _make_data((height, width), fill, backend, chunks)
     unit = "degree" if crs == 4326 else "m"
@@ -281,7 +291,7 @@ def from_template(name, resolution=None, *, preserve=None, backend="numpy",
         name=key,
         coords={"y": ys, "x": xs},
         dims=["y", "x"],
-        attrs={"res": (actual_res_x, actual_res_y), "crs": crs},
+        attrs={"res": (res_x, res_y), "crs": crs},
     )
     template["x"].attrs["units"] = unit
     template["y"].attrs["units"] = unit

xrspatial/tests/test_templates.py

@@ -82,6 +82,41 @@ def test_resolution_tuple():
     assert rx > ry
 
 
+def test_resolution_honored_exactly():
+    # the nyc bbox is 52814 m tall, not a whole number of 10 m cells, so res_y
+    # used to drift to ~10.0008. The far edge is nudged so res stays exact.
+    agg = from_template("nyc", resolution=10)
+    assert agg.attrs["res"] == (10.0, 10.0)
+
+
+def test_resolution_tuple_honored_exactly():
+    agg = from_template("conus", resolution=(10000, 5000))
+    assert agg.attrs["res"] == (10000.0, 5000.0)
+
+
+def test_coords_match_requested_resolution():
+    # pixel spacing equals the requested resolution on both axes
+    agg = from_template("nyc", resolution=10)
+    np.testing.assert_allclose(np.diff(agg.x.values), 10.0, atol=1e-6)
+    np.testing.assert_allclose(-np.diff(agg.y.values), 10.0, atol=1e-6)
+
+
+def test_nudge_keeps_centers_within_bbox():
+    # nudging the far edges out by < half a cell still leaves every pixel
+    # center inside the registry bbox.
+    agg = from_template("nyc", resolution=10)
+    left, bottom, right, top = _REGIONS["nyc"]["bounds"]
+    assert left <= agg.x.values.min() and agg.x.values.max() <= right
+    assert bottom <= agg.y.values.min() and agg.y.values.max() <= top
+
+
+def test_country_resolution_honored_exactly():
+    # country codes come back in EPSG:4326 (degrees) but go through the same
+    # nudge math, so a degree resolution is honored exactly too.
+    agg = from_template("FRA", resolution=0.25)
+    assert agg.attrs["res"] == (0.25, 0.25)
+
+
 def test_fill_and_dtype():
     agg = from_template("world")
     assert agg.dtype == np.float32
@@ -144,6 +179,10 @@ def test_dask_numpy_backend():
     np.testing.assert_array_equal(agg.x.values, ref.x.values)
     np.testing.assert_array_equal(agg.y.values, ref.y.values)
     assert agg.attrs == ref.attrs
+    # resolution is honored exactly on the dask path too
+    assert from_template("nyc", resolution=10, backend="dask+numpy").attrs[
+        "res"
+    ] == (10.0, 10.0)
     # values match once computed
     assert np.isnan(agg.compute().values).all()
 
@@ -264,6 +303,11 @@ def test_preserve_resolution_control():
     np.testing.assert_allclose(coarse.attrs["res"][0], 50000, rtol=0.05)
 
 
+def test_preserve_resolution_honored_exactly():
+    agg = from_template("conus", preserve="area", resolution=50000)
+    assert agg.attrs["res"] == (50000.0, 50000.0)
+
+
 def test_preserve_none_unchanged():
     a = from_template("conus")
     b = from_template("conus", preserve=None)