diff --git a/docs/src/whatsnew/latest.rst b/docs/src/whatsnew/latest.rst index 3493c5fec0..6675636685 100644 --- a/docs/src/whatsnew/latest.rst +++ b/docs/src/whatsnew/latest.rst @@ -36,8 +36,8 @@ This document explains the changes made to Iris for this release #. `@bjlittle`_, `@pp-mo`_ and `@trexfeathers`_ added support for unstructured meshes, as described by `UGRID`_. This involved adding a data model (:pull:`3968`, :pull:`4014`, :pull:`4027`, :pull:`4036`, :pull:`4053`) and API (:pull:`4063`, - :pull:`4064`), and supporting representation (:pull:`4033`, :pull:`4054`) and - loading (:pull:`4058`) of data on meshes. + :pull:`4064`), and supporting representation (:pull:`4033`, :pull:`4054`) of + data on meshes. Most of this new API can be found in :mod:`iris.experimental.ugrid`. The key objects introduced are :class:`iris.experimental.ugrid.Mesh`, :class:`iris.experimental.ugrid.MeshCoord` and @@ -49,9 +49,21 @@ This document explains the changes made to Iris for this release property :attr:`~iris.cube.Cube.mesh` which returns a :class:`~iris.experimental.ugrid.Mesh` if one is attached to the :class:`~iris.cube.Cube` via a :class:`~iris.experimental.ugrid.MeshCoord`. - Finally, the context manager :obj:`~iris.experimental.ugrid.PARSE_UGRID_ON_LOAD` + +#. `@trexfeathers`_ added support for loading unstructured mesh data from netcdf data, + for files using the `UGRID`_ conventions. + The context manager :obj:`~iris.experimental.ugrid.PARSE_UGRID_ON_LOAD` provides a way to load UGRID files so that :class:`~iris.cube.Cube`\ s can be returned with a :class:`~iris.experimental.ugrid.Mesh` attached. + (:pull:`4058`). + +#. `@pp-mo`_ added support to save cubes with meshes to netcdf files, using the + `UGRID`_ conventions. + The existing :meth:`iris.save` function now does this, when saving cubes with meshes. + A routine :meth:`iris.experimental.ugrid.save_mesh` allows saving + :class:`~iris.experimental.ugrid.Mesh` objects to netcdf *without* any associated data + (i.e. not attached to cubes). + (:pull:`4318` and :pull:`4339`). 🐛 Bugs Fixed diff --git a/lib/iris/experimental/ugrid/__init__.py b/lib/iris/experimental/ugrid/__init__.py index e31de8dd59..29fb1db6fa 100644 --- a/lib/iris/experimental/ugrid/__init__.py +++ b/lib/iris/experimental/ugrid/__init__.py @@ -45,11 +45,11 @@ from ...util import guess_coord_axis __all__ = [ - "CFUGridReader", "Connectivity", "ConnectivityMetadata", "load_mesh", "load_meshes", + "save_mesh", "Mesh", "Mesh1DConnectivities", "Mesh1DCoords", @@ -3729,9 +3729,8 @@ def _build_aux_coord(coord_var, file_path): Construct a :class:`~iris.coords.AuxCoord` from a given :class:`CFUGridAuxiliaryCoordinateVariable`, and guess its mesh axis. - todo: integrate with standard loading API post-pyke. - """ + # TODO: integrate with standard saving API when no longer 'experimental'. assert isinstance(coord_var, CFUGridAuxiliaryCoordinateVariable) attributes = {} attr_units = get_attr_units(coord_var, attributes) @@ -3745,7 +3744,8 @@ def _build_aux_coord(coord_var, file_path): # Fetch climatological - not allowed for a Mesh, but loading it will # mean an informative error gets raised. climatological = False - # TODO: use CF_ATTR_CLIMATOLOGY once re-integrated post-pyke. + # TODO: use CF_ATTR_CLIMATOLOGY on re-integration, when no longer + # 'experimental'. attr_climatology = getattr(coord_var, "climatology", None) if attr_climatology is not None: climatology_vars = coord_var.cf_group.climatology @@ -3781,9 +3781,8 @@ def _build_connectivity(connectivity_var, file_path, location_dims): :class:`CFUGridConnectivityVariable`, and identify the name of its first dimension. - todo: integrate with standard loading API post-pyke. - """ + # TODO: integrate with standard saving API when no longer 'experimental'. assert isinstance(connectivity_var, CFUGridConnectivityVariable) attributes = {} attr_units = get_attr_units(connectivity_var, attributes) @@ -3823,9 +3822,8 @@ def _build_mesh(cf, mesh_var, file_path): """ Construct a :class:`Mesh` from a given :class:`CFUGridMeshVariable`. - todo: integrate with standard loading API post-pyke. - """ + # TODO: integrate with standard saving API when no longer 'experimental'. assert isinstance(mesh_var, CFUGridMeshVariable) attributes = {} attr_units = get_attr_units(mesh_var, attributes) @@ -3966,9 +3964,8 @@ def _build_mesh_coords(mesh, cf_var): Construct a tuple of :class:`MeshCoord` using from a given :class:`Mesh` and :class:`~iris.fileformats.cf.CFVariable`. - todo: integrate with standard loading API post-pyke. - """ + # TODO: integrate with standard saving API when no longer 'experimental'. # Identify the cube's mesh dimension, for attaching MeshCoords. locations_dimensions = { "node": mesh.node_dimension, @@ -3985,3 +3982,55 @@ def _build_mesh_coords(mesh, cf_var): # END of loading section. ############################################################################### + + +############################################################################### +# SAVING + + +def save_mesh(mesh, filename, netcdf_format="NETCDF4"): + """ + Save mesh(es) to a netCDF file. + + Args: + + * mesh (:class:`iris.experimental.ugrid.Mesh` or iterable): + mesh(es) to save. + + * filename (string): + Name of the netCDF file to create. + + Kwargs: + + * netcdf_format (string): + Underlying netCDF file format, one of 'NETCDF4', 'NETCDF4_CLASSIC', + 'NETCDF3_CLASSIC' or 'NETCDF3_64BIT'. Default is 'NETCDF4' format. + + """ + # TODO: integrate with standard saving API when no longer 'experimental'. + + if isinstance(mesh, Iterable): + meshes = mesh + else: + meshes = [mesh] + + # Initialise Manager for saving + with netcdf.Saver(filename, netcdf_format) as sman: + # Iterate through the list. + for mesh in meshes: + # Get suitable dimension names. + mesh_dimensions, _ = sman._get_dim_names(mesh) + + # Create dimensions. + sman._create_cf_dimensions( + cube=None, dimension_names=mesh_dimensions + ) + + # Create the mesh components. + sman._add_mesh(mesh) + + # Add a conventions attribute. + # TODO: add 'UGRID' to conventions, when this is agreed with CF ? + sman.update_global_attributes( + Conventions=netcdf.CF_CONVENTIONS_VERSION + ) diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py index 6aed38af68..659632d7b4 100644 --- a/lib/iris/fileformats/netcdf.py +++ b/lib/iris/fileformats/netcdf.py @@ -13,6 +13,7 @@ """ import collections +import collections.abc from itertools import repeat, zip_longest import os import os.path @@ -1356,16 +1357,17 @@ def _create_cf_dimensions( """ unlimited_dim_names = [] - for coord in unlimited_dimensions: - try: - coord = cube.coord(name_or_coord=coord, dim_coords=True) - except iris.exceptions.CoordinateNotFoundError: - # coordinate isn't used for this cube, but it might be - # used for a different one - pass - else: - dim_name = self._get_coord_variable_name(cube, coord) - unlimited_dim_names.append(dim_name) + if unlimited_dimensions is not None: + for coord in unlimited_dimensions: + try: + coord = cube.coord(name_or_coord=coord, dim_coords=True) + except iris.exceptions.CoordinateNotFoundError: + # coordinate isn't used for this cube, but it might be + # used for a different one + pass + else: + dim_name = self._get_coord_variable_name(cube, coord) + unlimited_dim_names.append(dim_name) for dim_name in dimension_names: if dim_name not in self._dataset.dimensions: @@ -1375,7 +1377,7 @@ def _create_cf_dimensions( size = self._existing_dim[dim_name] self._dataset.createDimension(dim_name, size) - def _add_mesh(self, cube): + def _add_mesh(self, cube_or_mesh): """ Add the cube's mesh, and all related variables to the dataset. Includes all the mesh-element coordinate and connectivity variables. @@ -1387,8 +1389,9 @@ def _add_mesh(self, cube): Args: - * cube (:class:`iris.cube.Cube`): - A :class:`iris.cube.Cube` to be saved to a netCDF file. + * cube_or_mesh (:class:`iris.cube.Cube` + or :class:`iris.experimental.ugrid.Mesh`): + The Cube or Mesh being saved to the netCDF file. Returns: * cf_mesh_name (string or None): @@ -1397,7 +1400,17 @@ def _add_mesh(self, cube): """ cf_mesh_name = None - mesh = cube.mesh + + # Do cube- or -mesh-based save + from iris.cube import Cube + + if isinstance(cube_or_mesh, Cube): + cube = cube_or_mesh + mesh = cube.mesh + else: + cube = None # The underlying routines must support this ! + mesh = cube_or_mesh + if mesh: cf_mesh_name = self._name_coord_map.name(mesh) if cf_mesh_name is None: @@ -1421,7 +1434,7 @@ def _add_mesh(self, cube): if coord is None: continue # an awkward thing that mesh.coords does coord_name = self._create_generic_cf_array_var( - cube, + cube_or_mesh, [], coord, element_dims=(mesh_dims[location],), @@ -1459,7 +1472,11 @@ def _add_mesh(self, cube): # Has the 'other' dimension order, =reversed conn_dims = conn_dims[::-1] cf_conn_name = self._create_generic_cf_array_var( - cube, [], conn, element_dims=conn_dims, fill_value=-1 + cube_or_mesh, + [], + conn, + element_dims=conn_dims, + fill_value=-1, ) # Add essential attributes to the Connectivity variable. cf_conn_var = self._dataset.variables[cf_conn_name] @@ -1710,14 +1727,15 @@ def _add_aux_factories(self, cube, cf_var_cube, dimension_names): _setncattr(cf_var, "axis", "Z") _setncattr(cf_var, "formula_terms", formula_terms) - def _get_dim_names(self, cube): + def _get_dim_names(self, cube_or_mesh): """ Determine suitable CF-netCDF data dimension names. Args: - * cube (:class:`iris.cube.Cube`): - A :class:`iris.cube.Cube` to be saved to a netCDF file. + * cube_or_mesh (:class:`iris.cube.Cube` + or :class:`iris.experimental.ugrid.Mesh`): + The Cube or Mesh being saved to the netCDF file. Returns: mesh_dimensions, cube_dimensions @@ -1766,9 +1784,20 @@ def record_dimension(names_list, dim_name, length, matching_coords=[]): # Add the latest name to the list passed in. names_list.append(dim_name) + # Choose cube or mesh saving. + from iris.cube import Cube + + if isinstance(cube_or_mesh, Cube): + # there is a mesh, only if the cube has one + cube = cube_or_mesh + mesh = cube.mesh + else: + # there is no cube, only a mesh + cube = None + mesh = cube_or_mesh + # Get info on mesh, first. mesh_dimensions = [] - mesh = cube.mesh if mesh is None: cube_mesh_dim = None else: @@ -1835,58 +1864,65 @@ def record_dimension(names_list, dim_name, length, matching_coords=[]): # Store the mesh dims indexed by location mesh_location_dimnames[location] = dim_name - # Finally, identify the cube dimension which maps to the mesh: this - # is used below to recognise the mesh among the cube dimensions. - any_mesh_coord = cube.coords(mesh_coords=True)[0] - (cube_mesh_dim,) = any_mesh_coord.cube_dims(cube) + if cube is None: + cube_mesh_dim = None + else: + # Finally, identify the cube dimension which maps to the mesh: this + # is used below to recognise the mesh among the cube dimensions. + any_mesh_coord = cube.coords(mesh_coords=True)[0] + (cube_mesh_dim,) = any_mesh_coord.cube_dims(cube) + # Record actual file dimension names for each mesh saved. self._mesh_dims[mesh] = mesh_location_dimnames # Get the cube dimensions, in order. cube_dimensions = [] - for dim in range(cube.ndim): - if dim == cube_mesh_dim: - # Handle a mesh dimension: we already named this. - dim_coords = [] - dim_name = self._mesh_dims[mesh][cube.location] - else: - # Get a name from the dim-coord (if any). - dim_coords = cube.coords(dimensions=dim, dim_coords=True) - if dim_coords: - # Derive a dim name from a coord. - coord = dim_coords[0] # always have at least one - - # Add only dimensions that have not already been added. - dim_name = self._dim_names_and_coords.name(coord) - if dim_name is None: - # Not already present : create a unique dimension name - # from the coord. - dim_name = self._get_coord_variable_name(cube, coord) - while ( - dim_name in self._existing_dim - or dim_name in self._name_coord_map.names - ): - dim_name = self._increment_name(dim_name) - + if cube is not None: + for dim in range(cube.ndim): + if dim == cube_mesh_dim: + # Handle a mesh dimension: we already named this. + dim_coords = [] + dim_name = self._mesh_dims[mesh][cube.location] else: - # No CF-netCDF coordinates describe this data dimension. - # Make up a new, distinct dimension name - dim_name = f"dim{dim}" - if dim_name in self._existing_dim: - # Increment name if conflicted with one already existing. - if self._existing_dim[dim_name] != cube.shape[dim]: + # Get a name from the dim-coord (if any). + dim_coords = cube.coords(dimensions=dim, dim_coords=True) + if dim_coords: + # Derive a dim name from a coord. + coord = dim_coords[0] # always have at least one + + # Add only dimensions that have not already been added. + dim_name = self._dim_names_and_coords.name(coord) + if dim_name is None: + # Not already present : create a unique dimension name + # from the coord. + dim_name = self._get_coord_variable_name( + cube, coord + ) while ( dim_name in self._existing_dim - and self._existing_dim[dim_name] - != cube.shape[dim] or dim_name in self._name_coord_map.names ): dim_name = self._increment_name(dim_name) - # Record the dimension. - record_dimension( - cube_dimensions, dim_name, cube.shape[dim], dim_coords - ) + else: + # No CF-netCDF coordinates describe this data dimension. + # Make up a new, distinct dimension name + dim_name = f"dim{dim}" + if dim_name in self._existing_dim: + # Increment name if conflicted with one already existing. + if self._existing_dim[dim_name] != cube.shape[dim]: + while ( + dim_name in self._existing_dim + and self._existing_dim[dim_name] + != cube.shape[dim] + or dim_name in self._name_coord_map.names + ): + dim_name = self._increment_name(dim_name) + + # Record the dimension. + record_dimension( + cube_dimensions, dim_name, cube.shape[dim], dim_coords + ) return mesh_dimensions, cube_dimensions @@ -2042,14 +2078,15 @@ def _get_cube_variable_name(self, cube): cf_name = self.cf_valid_var_name(cf_name) return cf_name - def _get_coord_variable_name(self, cube, coord): + def _get_coord_variable_name(self, cube_or_mesh, coord): """ Returns a CF-netCDF variable name for a given coordinate-like element. Args: - * cube (:class:`iris.cube.Cube`): - The cube that contains the given coordinate. + * cube_or_mesh (:class:`iris.cube.Cube` + or :class:`iris.experimental.ugrid.Mesh`): + The Cube or Mesh being saved to the netCDF file. * coord (:class:`iris.coords._DimensionalMetadata`): An instance of a coordinate (or similar), for which a CF-netCDF variable name is required. @@ -2058,13 +2095,23 @@ def _get_coord_variable_name(self, cube, coord): A CF-netCDF variable name as a string. """ + # Support cube or mesh save. + from iris.cube import Cube + + if isinstance(cube_or_mesh, Cube): + cube = cube_or_mesh + mesh = cube.mesh + else: + cube = None + mesh = cube_or_mesh + if coord.var_name is not None: cf_name = coord.var_name else: name = coord.standard_name or coord.long_name if not name or set(name).intersection(string.whitespace): # We need to invent a name, based on its associated dimensions. - if cube.coords(coord): + if cube is not None and cube.coords(coord): # It is a regular cube coordinate. # Auto-generate a name based on the dims. name = "" @@ -2084,7 +2131,7 @@ def _get_coord_variable_name(self, cube, coord): assert isinstance(coord, Connectivity) location = coord.cf_role.split("_")[0] location_dim_attr = f"{location}_dimension" - name = getattr(cube.mesh, location_dim_attr) + name = getattr(mesh, location_dim_attr) # Convert to lower case and replace whitespace by underscores. cf_name = "_".join(name.lower().split()) @@ -2193,7 +2240,12 @@ def _set_cf_var_attributes(self, cf_var, element): _setncattr(cf_var, name, value) def _create_generic_cf_array_var( - self, cube, cube_dim_names, element, element_dims=None, fill_value=None + self, + cube_or_mesh, + cube_dim_names, + element, + element_dims=None, + fill_value=None, ): """ Create the associated CF-netCDF variable in the netCDF dataset for the @@ -2205,8 +2257,9 @@ def _create_generic_cf_array_var( Args: - * cube (:class:`iris.cube.Cube`): - The associated cube being saved to CF-netCDF file. + * cube_or_mesh (:class:`iris.cube.Cube` + or :class:`iris.experimental.ugrid.Mesh`): + The Cube or Mesh being saved to the netCDF file. * cube_dim_names (list of string): The name of each dimension of the cube. * element: @@ -2230,8 +2283,17 @@ def _create_generic_cf_array_var( The name of the CF-netCDF variable created. """ + # Support cube or mesh save. + from iris.cube import Cube + + if isinstance(cube_or_mesh, Cube): + cube = cube_or_mesh + else: + cube = None + # Work out the var-name to use. - cf_name = self._get_coord_variable_name(cube, element) + # N.B. the only part of this routine that may use a mesh _or_ a cube. + cf_name = self._get_coord_variable_name(cube_or_mesh, element) while cf_name in self._dataset.variables: cf_name = self._increment_name(cf_name) @@ -2295,7 +2357,7 @@ def _create_generic_cf_array_var( fill_value = None # Check if this is a dim-coord. - is_dimcoord = element in cube.dim_coords + is_dimcoord = cube is not None and element in cube.dim_coords if isinstance(element, iris.coords.CellMeasure): # Disallow saving of *masked* cell measures. diff --git a/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py b/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py index 16b90a766c..dffe3cd56f 100644 --- a/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py +++ b/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py @@ -19,7 +19,7 @@ from iris import save from iris.coords import AuxCoord from iris.cube import Cube, CubeList -from iris.experimental.ugrid import Connectivity, Mesh +from iris.experimental.ugrid import Connectivity, Mesh, save_mesh from iris.tests.stock import realistic_4d XY_LOCS = ("x", "y") @@ -683,11 +683,6 @@ def check_save_mesh(self, mesh): Use a name unique to this testcase, to avoid any clashes. """ - # FOR NOW: build a cube around the mesh. - # TODO: test via a 'save_mesh' function, when available. - cube = make_cube(mesh, location="node") - # N.B. using 'node' as always present, unlike the default 'face' - # use 'result_path' to name the file after the test function tempfile_path = self.result_path(ext=".nc") # Create a file of that name, but discard the result path and put it @@ -695,7 +690,7 @@ def check_save_mesh(self, mesh): tempfile_path = self.temp_dir / Path(tempfile_path).name # Save data to the file. - save(cube, tempfile_path) + save_mesh(mesh, tempfile_path) return tempfile_path @@ -1173,6 +1168,96 @@ def test_connectivity_names(self): ) self.assertEqual(expected_names, result_names, fail_msg) + def _check_two_different_meshes(self, vars): + # there are exactly 2 meshes in the file + mesh_names = vars_meshnames(vars) + self.assertEqual(sorted(mesh_names), ["Mesh2d", "Mesh2d_0"]) + + # they use different dimensions + # mesh1 + self.assertEqual( + vars_meshdim(vars, "node", mesh_name="Mesh2d"), "Mesh2d_nodes" + ) + self.assertEqual( + vars_meshdim(vars, "face", mesh_name="Mesh2d"), "Mesh2d_faces" + ) + if "edge_coordinates" in vars["Mesh2d"]: + self.assertEqual( + vars_meshdim(vars, "edge", mesh_name="Mesh2d"), "Mesh2d_edge" + ) + + # mesh2 + self.assertEqual( + vars_meshdim(vars, "node", mesh_name="Mesh2d_0"), "Mesh2d_nodes_0" + ) + self.assertEqual( + vars_meshdim(vars, "face", mesh_name="Mesh2d_0"), "Mesh2d_faces_0" + ) + if "edge_coordinates" in vars["Mesh2d_0"]: + self.assertEqual( + vars_meshdim(vars, "edge", mesh_name="Mesh2d_0"), + "Mesh2d_edge_0", + ) + + # the relevant coords + connectivities are also distinct + # mesh1 + self.assertEqual(vars["node_x"][_VAR_DIMS], ["Mesh2d_nodes"]) + self.assertEqual(vars["face_x"][_VAR_DIMS], ["Mesh2d_faces"]) + self.assertEqual( + vars["mesh2d_faces"][_VAR_DIMS], + ["Mesh2d_faces", "Mesh2d_face_N_nodes"], + ) + if "edge_coordinates" in vars["Mesh2d"]: + self.assertEqual(vars["longitude"][_VAR_DIMS], ["Mesh2d_edge"]) + self.assertEqual( + vars["mesh2d_edge"][_VAR_DIMS], + ["Mesh2d_edge", "Mesh2d_edge_N_nodes"], + ) + + # mesh2 + self.assertEqual(vars["node_x_0"][_VAR_DIMS], ["Mesh2d_nodes_0"]) + self.assertEqual(vars["face_x_0"][_VAR_DIMS], ["Mesh2d_faces_0"]) + self.assertEqual( + vars["mesh2d_faces_0"][_VAR_DIMS], + ["Mesh2d_faces_0", "Mesh2d_0_face_N_nodes"], + ) + if "edge_coordinates" in vars["Mesh2d_0"]: + self.assertEqual(vars["longitude_0"][_VAR_DIMS], ["Mesh2d_edge_0"]) + self.assertEqual( + vars["mesh2d_edge_0"][_VAR_DIMS], + ["Mesh2d_edge_0", "Mesh2d_0_edge_N_nodes"], + ) + + def test_multiple_identical_meshes(self): + mesh1 = make_mesh() + mesh2 = make_mesh() + + # Save and snapshot the result + tempfile_path = self.check_save_mesh([mesh1, mesh2]) + dims, vars = scan_dataset(tempfile_path) + + # Check there are two independent meshes + self._check_two_different_meshes(vars) + + def test_multiple_different_meshes(self): + # Create 2 meshes with different faces, but same edges. + # N.B. they should *not* then share an edge dimension ! + mesh1 = make_mesh(n_faces=3, n_edges=2) + mesh2 = make_mesh(n_faces=4, n_edges=2) + + # Save and snapshot the result + tempfile_path = self.check_save_mesh([mesh1, mesh2]) + dims, vars = scan_dataset(tempfile_path) + + # Check there are two independent meshes + self._check_two_different_meshes(vars) + + # Check the dims are as expected + self.assertEqual(dims["Mesh2d_faces"], 3) + self.assertEqual(dims["Mesh2d_faces_0"], 4) + self.assertEqual(dims["Mesh2d_edge"], 2) + self.assertEqual(dims["Mesh2d_edge_0"], 2) + if __name__ == "__main__": tests.main()