diff --git a/lib/iris/fileformats/netcdf/saver.py b/lib/iris/fileformats/netcdf/saver.py index 4a2474ba9b..a06d4eb9b4 100644 --- a/lib/iris/fileformats/netcdf/saver.py +++ b/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. diff --git a/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py b/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py index df04016cb7..4985d819a3 100644 --- a/lib/iris/tests/integration/netcdf/derived_bounds/test_bounds_files.py +++ b/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 + ) diff --git a/lib/iris/tests/results/integration/netcdf/derived_bounds/TestBoundsFiles/test_save_primary_with_db.cdl b/lib/iris/tests/results/integration/netcdf/derived_bounds/TestBoundsFiles/test_save_primary_with_db.cdl new file mode 100644 index 0000000000..f5a0de6480 --- /dev/null +++ b/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" ; +} diff --git a/lib/iris/tests/results/integration/netcdf/derived_bounds/TestBoundsFiles/test_save_primary_without_db.cdl b/lib/iris/tests/results/integration/netcdf/derived_bounds/TestBoundsFiles/test_save_primary_without_db.cdl new file mode 100644 index 0000000000..f2879e22ea --- /dev/null +++ b/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" ; +}