from_template: use CF Conventions metadata instead of crs_units/crs_name

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,22 @@ 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. 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
     --------
@@ -302,14 +323,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 +334,13 @@ 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

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