properly decode and encode ellipsoids for the cf convention

docs/changelog.md

@@ -1,5 +1,10 @@
 # Changelog
 
+## 0.6.1 (_unreleased_)
+
+- Do not mutate the input dataset when decoding ({pull}`226`)
+- Properly decode ellipsoids for the CF convention ({pull}`228`)
+
 ## 0.6.0 (2026-02-05)
 
 - Support variable-sized cells with the `"zuniq"` indexing scheme in `healpix` ({pull}`207`)

pixi.toml

@@ -100,13 +100,14 @@ isort = "*"
 jupyterlab = "*"
 jupyter-resource-usage = "*"
 jupyterlab_code_formatter = "*"
-python-build = ">=1.3.0,<2"
-h3-py = ">=4.3.0,<5"
-geopandas = ">=1.1.1,<2"
-pyinstrument = ">=5.1.1,<6"
-jupyterlab-myst = ">=2.4.2,<3"
-h5netcdf = ">=1.7.3,<2"
-zarr = ">=3.1.5,<4"
+python-build = ">=1.3.0"
+h3-py = ">=4.3.0"
+geopandas = ">=1.1.1"
+pyinstrument = ">=5.1.1"
+jupyterlab-myst = ">=2.4.2"
+h5netcdf = ">=1.7.3"
+zarr = ">=3.1.5"
+cf-xarray = ">=0.11.1"
 
 [environments]
 nightly = { features = [

xdggs/conventions/cf.py

@@ -1,4 +1,4 @@
-from collections.abc import Hashable
+from collections.abc import Hashable, Sequence
 from typing import Any, Literal
 
 import numpy as np
@@ -22,6 +22,23 @@ def remove_grid_mapping(ds, name):
     return new
 
 
+def remove_keys(mapping: dict[str, Any], exclude: Sequence[str]) -> dict[str, Any]:
+    return {k: v for k, v in mapping.items() if k not in exclude}
+
+
+ellipsoid_attribute_translations = {
+    "semi_major_axis": "semimajor_axis",
+    "semi_minor_axis": "semiminor_axis",
+    "earth_radius": "radius",
+    "reference_ellipsoid_name": "name",
+    "inverse_flattening": "inverse_flattening",
+}
+
+
+def extract_ellipsoid_parameters(mapping: dict[str, Any]) -> dict[str, Any]:
+    return {k: v for k, v in mapping.items() if k in ellipsoid_attribute_translations}
+
+
 @register_convention("cf")
 class Cf(Convention):
     def translate_keys(
@@ -35,6 +52,17 @@ def translate_keys(
 
         return {translations.get(key, key): value for key, value in mapping.items()}
 
+    def translate_ellipsoid(
+        self,
+        mapping: dict[str, Any],
+        direction: Literal["forward", "inverse"] = "forward",
+    ) -> dict[str, Any] | None:
+        translations = ellipsoid_attribute_translations
+        if direction == "inverse":
+            translations = {v: k for k, v in translations.items()}
+
+        return {translations.get(key, key): value for key, value in mapping.items()}
+
     def decode(
         self,
         ds: xr.Dataset,
@@ -70,7 +98,13 @@ def decode(
                 )
             name = coord_names[0]
 
-        grid_info = self.translate_keys(crs.attrs, direction="forward")
+        ellipsoid = extract_ellipsoid_parameters(crs.attrs)
+        grid_info = self.translate_keys(
+            remove_keys(crs.attrs, ellipsoid), direction="forward"
+        )
+        if parsed_ellipsoid := self.translate_ellipsoid(ellipsoid, direction="forward"):
+            grid_info["ellipsoid"] = parsed_ellipsoid
+
         grid_name = grid_info["grid_name"]
 
         var = ds.variables[name].copy(deep=False)
@@ -100,8 +134,11 @@ def encode(self, ds: xr.Dataset, *, encoding: dict[str, Any] | None = None):
 
         grid_name = infer_grid_name(ds.dggs.index)
         grid_info_dict = grid_info.to_dict()
-        metadata = self.translate_keys(grid_info_dict, direction="inverse")
+        ellipsoid = grid_info_dict.pop("ellipsoid", {})
 
+        metadata = self.translate_keys(
+            grid_info_dict, direction="inverse"
+        ) | self.translate_ellipsoid(ellipsoid, direction="inverse")
         crs = xr.Variable((), np.int8(0), metadata)
 
         additional_var_attrs = {"grid_mapping": grid_name}

xdggs/conventions/xdggs.py

@@ -1,3 +1,4 @@
+import copy
 from collections.abc import Hashable
 from typing import Any
 
@@ -40,17 +41,24 @@ def decode(
         [dim] = var.dims
 
         if grid_info is None:
-            grid_info = var.attrs
+            grid_info = copy.deepcopy(var.attrs)
         elif isinstance(grid_info, DGGSInfo):
             # TODO: avoid serializing / deserializing cycle
             grid_info = grid_info.to_dict()
 
-        grid_name = grid_info["grid_name"]
+        try:
+            grid_name = grid_info["grid_name"]
+        except KeyError:
+            raise DecoderError(
+                "xdggs convention: no grid name found (`grid_name`)."
+                " Try choosing a different convention or verify the dataset metadata."
+            ) from None
+
         if grid_name not in GRID_REGISTRY:
             raise DecoderError(f"xdggs convention: unknown grid name: {grid_name}")
         index_cls = GRID_REGISTRY[grid_name]
 
-        var_ = var.copy(deep=True)
+        var_ = var.copy(deep=False)
         var_.attrs = grid_info
         index = index_cls.from_variables({name: var_}, options=index_options)
 
@@ -65,4 +73,4 @@ def encode(self, ds: xr.Dataset, *, encoding: None = None) -> xr.Dataset:
         # TODO: `assign_coords` + `assign_attrs` drops the index
         ds_ = ds.copy(deep=False)
         ds_[coord].attrs.update(metadata)
-        return ds_
+        return ds_.drop_indexes(coord)

xdggs/ellipsoid.py

@@ -29,7 +29,15 @@ class Ellipsoid:
 
     @classmethod
     def from_dict(cls, mapping):
-        return cls(**mapping)
+        semimajor = mapping["semimajor_axis"]
+        if (semiminor := mapping.get("semiminor_axis")) is not None:
+            inverse_flattening = semimajor / (semimajor - semiminor)
+        else:
+            inverse_flattening = mapping["inverse_flattening"]
+        name = mapping.get("name")
+        return cls(
+            semimajor_axis=semimajor, inverse_flattening=inverse_flattening, name=name
+        )
 
     def to_dict(self):
         mapping = asdict(self)

xdggs/h3.py

@@ -242,5 +242,8 @@ def grid_info(self) -> H3Info:
     def _replace(self, new_index: xr.Index):
         return type(self)(new_index, self._dim, self._name, self._grid)
 
+    def __repr__(self):
+        return f"<H3Index(level={self._grid.level})>"
+
     def _repr_inline_(self, max_width: int):
         return f"H3Index(level={self._grid.level})"

xdggs/healpix.py

@@ -190,6 +190,8 @@ def translate_nside(nside):
         def translate_ellipsoid(value):
             if isinstance(value, (str, Sphere, Ellipsoid)):
                 return value
+            elif value is None or not value:
+                return value
 
             return parse_ellipsoid(value)
 
@@ -735,5 +737,13 @@ def _replace(self, new_index: xr.Index):
     def grid_info(self) -> HealpixInfo:
         return self._grid
 
+    def __repr__(self):
+        return "\n".join(
+            [
+                f"<HealpixIndex(kind={self._kind})>",
+                repr(self._grid),
+            ]
+        )
+
     def _repr_inline_(self, max_width: int):
         return f"HealpixIndex(level={self._grid.level}, indexing_scheme={self._grid.indexing_scheme}, kind={self._kind})"