Save derived bounds

lib/iris/fileformats/netcdf/saver.py

@@ -29,6 +29,7 @@
 from dask.delayed import Delayed
 import numpy as np
 
+from iris import FUTURE
 from iris._deprecation import warn_deprecated
 from iris._lazy_data import _co_realise_lazy_arrays, is_lazy_data
 from iris.aux_factory import (
@@ -1068,11 +1069,13 @@ def _add_aux_factories(self, cube, cf_var_cube, dimension_names):
                 cf_name = self._name_coord_map.name(primary_coord)
                 cf_var = self._dataset.variables[cf_name]
 
-                names = {
-                    key: self._name_coord_map.name(coord)
-                    for key, coord in factory.dependencies.items()
+                term_varnames = {
+                    term: self._name_coord_map.name(coord)
+                    for term, coord in factory.dependencies.items()
                 }
-                formula_terms = factory_defn.formula_terms_format.format(**names)
+                formula_terms = factory_defn.formula_terms_format.format(
+                    **term_varnames
+                )
                 std_name = factory_defn.std_name
 
                 if hasattr(cf_var, "formula_terms"):
@@ -1114,6 +1117,39 @@ def _add_aux_factories(self, cube, cf_var_cube, dimension_names):
                     _setncattr(cf_var, "axis", "Z")
                     _setncattr(cf_var, "formula_terms", formula_terms)
 
+                if FUTURE.derived_bounds:
+                    # ensure that the primary variable *bounds*, if any, obey the CF
+                    #  encoding rule : the bounds variable of a parametric coordinate
+                    #  must itself have a "formula_terms" attribute.
+                    # See : https://cfconventions.org/Data/cf-conventions/cf-conventions-1.12/cf-conventions.html#boundaries-and-formula-terms
+                    bounds_varname = getattr(cf_var, "bounds", None)
+                    cf_bounds_var = self._dataset.variables.get(bounds_varname, None)
+                    if (
+                        cf_bounds_var is not None
+                        and getattr(cf_bounds_var, "formula_terms", None) is None
+                    ):
+                        # We need a bounds formula, and there is none already attached.
+                        # Construct and add one, mirroring the main formula.
+                        def boundsterm_varname(term_varname):
+                            """Identify the boundsvar for a given term var."""
+                            # First establish a fallback default, meaning 'no bounds'.
+                            result = term_varname
+                            # Follow links (if they exist) to find the bounds var.
+                            termvar = self._dataset.variables.get(term_varname)
+                            boundsname = getattr(termvar, "bounds", None)
+                            if boundsname in self._dataset.variables:
+                                result = boundsname
+                            return result
+
+                        boundsterm_varnames = {
+                            term: boundsterm_varname(term_varname)
+                            for term, term_varname in term_varnames.items()
+                        }
+                        bounds_formula_terms = factory_defn.formula_terms_format.format(
+                            **boundsterm_varnames
+                        )
+                        _setncattr(cf_bounds_var, "formula_terms", bounds_formula_terms)
+
     def _get_dim_names(self, cube_or_mesh):
         """Determine suitable CF-netCDF data dimension names.
 

lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py

@@ -18,13 +18,20 @@
 
 import iris
 from iris import FUTURE, sample_data_path
+from iris.tests._shared_utils import assert_CDL
 from iris.tests.stock.netcdf import ncgen_from_cdl
 
 db_testfile_path = Path(__file__).parent / "temp_nc_sources" / "a_new_file.nc"
 legacy_filepath = sample_data_path("hybrid_height.nc")
 
+_DERIVED_BOUNDS_IDS = ["with_db", "without_db"]
+_DERIVED_BOUNDS_VALUES = [True, False]
+_DERIVED_BOUNDS_ID_MAP = {
+    value: id for value, id in zip(_DERIVED_BOUNDS_VALUES, _DERIVED_BOUNDS_IDS)
+}
 
-@pytest.fixture(params=[True, False], ids=["with_db", "without_db"])
+
+@pytest.fixture(params=_DERIVED_BOUNDS_VALUES, ids=_DERIVED_BOUNDS_IDS)
 def derived_bounds(request):
     db = request.param
     with FUTURE.context(derived_bounds=db):
@@ -181,3 +188,35 @@ def test_load_primary_cf_style(derived_bounds, cf_primary_sample_path):
     assert co_P0.var_name == "P0"
     assert not co_P0.has_bounds()
     assert main_cube.coord_dims(co_P0) == ()
+
+
+@pytest.fixture()
+def tmp_ncdir(tmp_path_factory):
+    yield tmp_path_factory.mktemp("_temp_netcdf_dir")
+
+
+def test_save_primary_cf_style(
+    derived_bounds, cf_primary_sample_path, request, tmp_ncdir
+):
+    """Check how our 'standard primary encoded' derived coordinate example saves.
+
+    Test against saved snapshot CDL, with and without FUTURE.derived_bounds enabled.
+    """
+    # N.B. always **load** with derived bounds enabled, as the content implies it...
+    with FUTURE.context(derived_bounds=True):
+        test_cube = iris.load(cf_primary_sample_path, "air_temperature")
+
+    # ... but whether we **save** with full derived-bounds handling depends on test mode.
+    db_id = _DERIVED_BOUNDS_ID_MAP[derived_bounds]
+
+    nc_filename = f"test_save_primary_{db_id}.nc"
+    cdl_filename = nc_filename.replace("nc", "cdl")
+    nc_filepath = tmp_ncdir / nc_filename
+    cdl_filepath = "integration/netcdf/derived_bounds/TestBoundsFiles/" + cdl_filename
+
+    # Save to test netcdf file
+    iris.save(test_cube, nc_filepath)
+    # Dump to CDL, check against stored reference snapshot.
+    assert_CDL(
+        request=request, netcdf_filename=nc_filepath, reference_filename=cdl_filepath
+    )

lib/iris/tests/results/integration/netcdf/derived_bounds/TestBoundsFiles/test_save_primary_with_db.cdl

@@ -0,0 +1,44 @@
+dimensions:
+	bnds = 2 ;
+	dim1 = 1 ;
+	dim2 = 1 ;
+	eta = 1 ;
+variables:
+	float temp(eta, dim1, dim2) ;
+		temp:standard_name = "air_temperature" ;
+		temp:units = "K" ;
+		temp:coordinates = "A B P0 PS ap" ;
+	double eta(eta) ;
+		eta:axis = "Z" ;
+		eta:bounds = "eta_bnds" ;
+		eta:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
+		eta:long_name = "eta at full levels" ;
+		eta:positive = "down" ;
+	double eta_bnds(eta, bnds) ;
+	double P0 ;
+		P0:units = "Pa" ;
+	double PS(dim1, dim2) ;
+		PS:units = "Pa" ;
+	double A(eta) ;
+		A:bounds = "A_bnds" ;
+		A:units = "1" ;
+		A:long_name = "a coefficient for vertical coordinate at full levels" ;
+	double A_bnds(eta, bnds) ;
+	double B(eta) ;
+		B:bounds = "B_bnds" ;
+		B:units = "1" ;
+		B:long_name = "b coefficient for vertical coordinate at full levels" ;
+	double B_bnds(eta, bnds) ;
+	double ap(eta) ;
+		ap:bounds = "ap_bnds" ;
+		ap:units = "Pa" ;
+		ap:long_name = "vertical pressure" ;
+		ap:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
+		ap:axis = "Z" ;
+		ap:formula_terms = "ap: ap b: B ps: PS" ;
+	double ap_bnds(eta, bnds) ;
+		ap_bnds:formula_terms = "ap: ap_bnds b: B_bnds ps: PS" ;
+
+// global attributes:
+		:Conventions = "CF-1.7" ;
+}

lib/iris/tests/results/integration/netcdf/derived_bounds/TestBoundsFiles/test_save_primary_without_db.cdl

@@ -0,0 +1,43 @@
+dimensions:
+	bnds = 2 ;
+	dim1 = 1 ;
+	dim2 = 1 ;
+	eta = 1 ;
+variables:
+	float temp(eta, dim1, dim2) ;
+		temp:standard_name = "air_temperature" ;
+		temp:units = "K" ;
+		temp:coordinates = "A B P0 PS ap" ;
+	double eta(eta) ;
+		eta:axis = "Z" ;
+		eta:bounds = "eta_bnds" ;
+		eta:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
+		eta:long_name = "eta at full levels" ;
+		eta:positive = "down" ;
+	double eta_bnds(eta, bnds) ;
+	double P0 ;
+		P0:units = "Pa" ;
+	double PS(dim1, dim2) ;
+		PS:units = "Pa" ;
+	double A(eta) ;
+		A:bounds = "A_bnds" ;
+		A:units = "1" ;
+		A:long_name = "a coefficient for vertical coordinate at full levels" ;
+	double A_bnds(eta, bnds) ;
+	double B(eta) ;
+		B:bounds = "B_bnds" ;
+		B:units = "1" ;
+		B:long_name = "b coefficient for vertical coordinate at full levels" ;
+	double B_bnds(eta, bnds) ;
+	double ap(eta) ;
+		ap:bounds = "ap_bnds" ;
+		ap:units = "Pa" ;
+		ap:long_name = "vertical pressure" ;
+		ap:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
+		ap:axis = "Z" ;
+		ap:formula_terms = "ap: ap b: B ps: PS" ;
+	double ap_bnds(eta, bnds) ;
+
+// global attributes:
+		:Conventions = "CF-1.7" ;
+}