From 14353510168ca7c2332180ac664371b8221197b2 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Thu, 25 Jun 2026 20:03:05 -0700 Subject: [PATCH 1/2] Use CF Conventions metadata in from_template (#3531) Replace the ad-hoc crs_units/crs_name attrs with CF grid-mapping keys grid_mapping_name + crs_wkt (via pyproj.CRS.to_cf, best-effort) and move units onto the x/y coordinates with CF standard_names. EPSG:4326 axes use degrees_east/degrees_north; projected axes use metres. --- xrspatial/templates.py | 73 +++++++++++++++++++----------- xrspatial/tests/test_templates.py | 74 ++++++++++++++++++++----------- 2 files changed, 97 insertions(+), 50 deletions(-) diff --git a/xrspatial/templates.py b/xrspatial/templates.py index b18870bb9..0e03dd958 100644 --- a/xrspatial/templates.py +++ b/xrspatial/templates.py @@ -173,18 +173,29 @@ def _make_data(shape, fill, backend, chunks): ) -def _crs_name(crs): - """Best-effort human-readable name for an EPSG code. - - Resolves through the same two-tier path as the rest of the library - (built-in ``LiteCRS`` table first, pyproj fallback). Returns ``None`` - when the code is outside the built-in table and pyproj is not - installed, so the default (non-reproject) path stays dependency-free. +def _cf_crs_attrs(crs): + """CF Conventions grid-mapping attributes for an EPSG code. + + Emits the two canonical CF-1.7 identifiers via :meth:`pyproj.CRS.to_cf`: + ``grid_mapping_name`` (the controlled-vocabulary projection token, e.g. + ``'albers_conical_equal_area'``) and ``crs_wkt`` (the full WKT, which + carries the human-readable CRS name). This mirrors the rest of the + library, which identifies a CRS by ``crs`` / ``crs_wkt`` and derives + anything else with pyproj on demand. + + Best-effort: without pyproj installed the mapping is empty, so the + default (non-reproject) path stays dependency-free, same as the old + ``crs_name`` behaviour. ``grid_mapping_name`` is also left out for + projections CF does not define (e.g. Equal Earth), where ``to_cf`` + returns ``crs_wkt`` alone. """ try: - return _resolve_crs(crs).name - except (ImportError, ValueError, TypeError): - return None + import pyproj + except ImportError: + return {} + cf = pyproj.CRS.from_epsg(int(crs)).to_cf() + return {k: cf[k] for k in ("grid_mapping_name", "crs_wkt") + if cf.get(k) is not None} def from_template(name, resolution=None, *, preserve=None, backend="numpy", @@ -233,12 +244,20 @@ def from_template(name, resolution=None, *, preserve=None, backend="numpy", Returns ------- template : xarray.DataArray - Empty 2-D raster with ``dims=('y', 'x')``, pixel-center coordinates, - and ``attrs`` carrying ``res``, ``crs``, ``crs_units`` (``'m'`` or - ``'degree'``), and ``crs_name`` (the projection's human-readable - name, e.g. ``'NAD83 / Conus Albers'``). ``crs_name`` is omitted only - when the EPSG code is outside the built-in table and pyproj is not - installed. + Empty 2-D raster with ``dims=('y', 'x')`` and pixel-center + coordinates. ``attrs`` carries ``res`` and ``crs`` plus the CF + Conventions grid-mapping keys ``grid_mapping_name`` (the projection + token, e.g. ``'albers_conical_equal_area'``) and ``crs_wkt`` (full + WKT, which carries the human-readable CRS name). The ``x`` / ``y`` + coordinates follow CF axis conventions: ``units`` ``'m'`` with + ``standard_name`` ``'projection_x_coordinate'`` / + ``'projection_y_coordinate'`` for projected templates, or + ``'degrees_east'`` / ``'degrees_north'`` with ``standard_name`` + ``'longitude'`` / ``'latitude'`` for EPSG:4326. The grid-mapping keys + require pyproj; without it they are omitted (the default, + dependency-free path), and ``grid_mapping_name`` is also omitted for + projections CF does not define (e.g. Equal Earth), leaving + ``crs_wkt`` alone. Examples -------- @@ -302,14 +321,9 @@ def from_template(name, resolution=None, *, preserve=None, backend="numpy", ys, xs = _make_output_coords((left, bottom, right, top), (height, width)) data = _make_data((height, width), fill, backend, chunks) - # Every template emits either EPSG:4326 (country codes, degrees) or a - # metre-based projected CRS (curated regions and all preserve paths). - unit = "degree" if crs == 4326 else "m" - attrs = {"res": (res_x, res_y), "crs": crs, "crs_units": unit} - crs_name = _crs_name(crs) - if crs_name is not None: - attrs["crs_name"] = crs_name + attrs = {"res": (res_x, res_y), "crs": crs} + attrs.update(_cf_crs_attrs(crs)) template = xr.DataArray( data, @@ -318,6 +332,15 @@ def from_template(name, resolution=None, *, preserve=None, backend="numpy", dims=["y", "x"], attrs=attrs, ) - template["x"].attrs["units"] = unit - template["y"].attrs["units"] = unit + # CF coordinate metadata (CF Conventions sec. 4). Every template emits + # either EPSG:4326 (country codes) or a metre-based projected CRS + # (curated regions and all preserve paths). + if crs == 4326: + template["x"].attrs.update(units="degrees_east", standard_name="longitude") + template["y"].attrs.update(units="degrees_north", standard_name="latitude") + else: + template["x"].attrs.update( + units="m", standard_name="projection_x_coordinate") + template["y"].attrs.update( + units="m", standard_name="projection_y_coordinate") return template diff --git a/xrspatial/tests/test_templates.py b/xrspatial/tests/test_templates.py index 3a17e3d7a..767852edf 100644 --- a/xrspatial/tests/test_templates.py +++ b/xrspatial/tests/test_templates.py @@ -55,7 +55,7 @@ def test_nyc_resolves(): def test_country_code(): agg = from_template("FRA") assert agg.attrs["crs"] == 4326 - assert agg.x.attrs["units"] == "degree" + assert agg.x.attrs["units"] == "degrees_east" assert np.isnan(agg.values).all() assert agg.shape[0] > 1 and agg.shape[1] > 1 assert agg.name == "FRA" @@ -341,43 +341,67 @@ def test_preserve_downstream_slope(): assert out.shape == agg.shape -def test_crs_units_attr(): - # crs_units mirrors the coordinate units: metres for a projected - # template, degrees for a country code in EPSG:4326. +def test_cf_coordinate_units(): + # Units live on the coordinates (CF Conventions sec. 4), not on a + # crs_units attr. Projected templates use metres; EPSG:4326 uses the + # CF degree spellings. proj = from_template("conus") - assert proj.attrs["crs_units"] == "m" == proj.x.attrs["units"] + assert "crs_units" not in proj.attrs + assert proj.x.attrs["units"] == "m" + assert proj.x.attrs["standard_name"] == "projection_x_coordinate" + assert proj.y.attrs["units"] == "m" + assert proj.y.attrs["standard_name"] == "projection_y_coordinate" + geo = from_template("FRA") - assert geo.attrs["crs_units"] == "degree" == geo.x.attrs["units"] + assert "crs_units" not in geo.attrs + assert geo.x.attrs["units"] == "degrees_east" + assert geo.x.attrs["standard_name"] == "longitude" + assert geo.y.attrs["units"] == "degrees_north" + assert geo.y.attrs["standard_name"] == "latitude" -def test_crs_name_attr(): - # 5070 is in the built-in LiteCRS table, so the name is deterministic - # regardless of whether pyproj is installed. - assert from_template("conus").attrs["crs_name"] == "NAD83 / Conus Albers" - assert from_template("FRA").attrs["crs_name"] == "WGS 84" +def test_cf_grid_mapping_attrs(): + # crs_name is gone; the CF grid-mapping keys identify the projection. + proj = from_template("conus") + assert "crs_name" not in proj.attrs + assert proj.attrs["grid_mapping_name"] == "albers_conical_equal_area" + assert "NAD83 / Conus Albers" in proj.attrs["crs_wkt"] + geo = from_template("FRA") + assert geo.attrs["grid_mapping_name"] == "latitude_longitude" + assert "WGS 84" in geo.attrs["crs_wkt"] -def test_crs_name_preserve_path(): - # The preserve path requires pyproj, so the chosen projection's name - # is always available. - name = from_template("conus", preserve="area").attrs["crs_name"] - assert name == "NAD83 / Conus Albers" +def test_cf_grid_mapping_preserve_path(): + # The preserve path requires pyproj, so the CF keys are always present. + agg = from_template("conus", preserve="area") + assert agg.attrs["grid_mapping_name"] == "albers_conical_equal_area" + assert "Conus Albers" in agg.attrs["crs_wkt"] -def test_crs_name_omitted_without_pyproj(monkeypatch): - # When the EPSG code is outside the built-in table and pyproj can't be - # imported, crs_name is left off rather than raising. - import xrspatial.templates as templates - def _no_crs(_crs): - raise ImportError("pyproj not installed") +def test_grid_mapping_omitted_for_equal_earth(): + # Equal Earth (the preserve='area' fallback for the world bbox) has no + # CF grid mapping, so grid_mapping_name is left off and crs_wkt stands + # alone. + agg = from_template("world", preserve="area") + assert agg.attrs["crs"] == 8857 + assert "grid_mapping_name" not in agg.attrs + assert "Equal Earth" in agg.attrs["crs_wkt"] - monkeypatch.setattr(templates, "_resolve_crs", _no_crs) + +def test_cf_attrs_omitted_without_pyproj(monkeypatch): + # Without pyproj the default (non-reproject) path stays dependency-free: + # the CF grid-mapping keys are left off rather than raising. + import sys + + monkeypatch.setitem(sys.modules, "pyproj", None) agg = from_template("conus") - assert "crs_name" not in agg.attrs + assert "grid_mapping_name" not in agg.attrs + assert "crs_wkt" not in agg.attrs # The rest of the contract still holds. assert agg.attrs["crs"] == 5070 - assert agg.attrs["crs_units"] == "m" + assert agg.x.attrs["units"] == "m" + assert agg.x.attrs["standard_name"] == "projection_x_coordinate" def test_lite_crs_name_property(): From 24ea99fabc723e862da64bd03a5e3de0d3c61007 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Thu, 25 Jun 2026 20:06:07 -0700 Subject: [PATCH 2/2] Address review nits: consistent wrapping + CF flat-attrs docstring caveat (#3531) --- xrspatial/templates.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/xrspatial/templates.py b/xrspatial/templates.py index 0e03dd958..a333e0969 100644 --- a/xrspatial/templates.py +++ b/xrspatial/templates.py @@ -257,7 +257,9 @@ def from_template(name, resolution=None, *, preserve=None, backend="numpy", require pyproj; without it they are omitted (the default, dependency-free path), and ``grid_mapping_name`` is also omitted for projections CF does not define (e.g. Equal Earth), leaving - ``crs_wkt`` alone. + ``crs_wkt`` alone. These keys sit directly on ``attrs`` (the library's + flat CRS convention), not on a separate CF grid-mapping variable, so a + strict CF reader will not auto-detect them as a grid mapping. Examples -------- @@ -339,8 +341,6 @@ def from_template(name, resolution=None, *, preserve=None, backend="numpy", template["x"].attrs.update(units="degrees_east", standard_name="longitude") template["y"].attrs.update(units="degrees_north", standard_name="latitude") else: - template["x"].attrs.update( - units="m", standard_name="projection_x_coordinate") - template["y"].attrs.update( - units="m", standard_name="projection_y_coordinate") + template["x"].attrs.update(units="m", standard_name="projection_x_coordinate") + template["y"].attrs.update(units="m", standard_name="projection_y_coordinate") return template