diff --git a/xrspatial/templates.py b/xrspatial/templates.py index 7e01efa5e..c82a46f62 100644 --- a/xrspatial/templates.py +++ b/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 diff --git a/xrspatial/tests/test_templates.py b/xrspatial/tests/test_templates.py index 37d3d6f6b..99ba0fe7c 100644 --- a/xrspatial/tests/test_templates.py +++ b/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)