From 957e84fb6b74d750d176827c4aea4bc34ccb9fd0 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Fri, 27 Aug 2021 09:37:40 +0100 Subject: [PATCH 01/26] Refactor onto a common multidim var routine. --- lib/iris/fileformats/netcdf.py | 295 ++++++------------ .../unit/fileformats/netcdf/test_Saver.py | 2 +- 2 files changed, 104 insertions(+), 193 deletions(-) diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py index af1c764247..4913ffe969 100644 --- a/lib/iris/fileformats/netcdf.py +++ b/lib/iris/fileformats/netcdf.py @@ -1361,7 +1361,6 @@ def _add_inner_related_vars( cf_var_cube, dimension_names, coordlike_elements, - saver_create_method, role_attribute_name, ): # Common method to create a set of file variables and attach them to @@ -1373,7 +1372,9 @@ def _add_inner_related_vars( ): # Create the associated CF-netCDF variable. if element not in self._name_coord_map.coords: - cf_name = saver_create_method(cube, dimension_names, element) + cf_name = self._create_generic_cf_array_var( + cube, dimension_names, element + ) self._name_coord_map.append(cf_name, element) else: cf_name = self._name_coord_map.name(element) @@ -1409,7 +1410,6 @@ def _add_aux_coords(self, cube, cf_var_cube, dimension_names): cf_var_cube, dimension_names, cube.aux_coords, - self._create_cf_coord_variable, "coordinates", ) @@ -1432,7 +1432,6 @@ def _add_cell_measures(self, cube, cf_var_cube, dimension_names): cf_var_cube, dimension_names, cube.cell_measures(), - self._create_cf_cell_measure_variable, "cell_measures", ) @@ -1456,7 +1455,6 @@ def _add_ancillary_variables(self, cube, cf_var_cube, dimension_names): cf_var_cube, dimension_names, cube.ancillary_variables(), - self._create_cf_ancildata_variable, "ancillary_variables", ) @@ -1476,7 +1474,7 @@ def _add_dim_coords(self, cube, dimension_names): for coord in cube.dim_coords: # Create the associated coordinate CF-netCDF variable. if coord not in self._name_coord_map.coords: - cf_name = self._create_cf_coord_variable( + cf_name = self._create_generic_cf_array_var( cube, dimension_names, coord ) self._name_coord_map.append(cf_name, coord) @@ -1554,7 +1552,7 @@ def _add_aux_factories(self, cube, cf_var_cube, dimension_names): name = self._formula_terms_cache.get(key) if name is None: # Create a new variable - name = self._create_cf_coord_variable( + name = self._create_generic_cf_array_var( cube, dimension_names, primary_coord ) cf_var = self._dataset.variables[name] @@ -1736,7 +1734,7 @@ def _create_cf_bounds(self, coord, cf_var, cf_name): None """ - if coord.has_bounds(): + if hasattr(coord, "has_bounds") and coord.has_bounds(): # Get the values in a form which is valid for the file format. bounds = self._ensure_valid_dtype( coord.bounds, "the bounds of coordinate", coord @@ -1793,15 +1791,15 @@ def _get_cube_variable_name(self, cube): def _get_coord_variable_name(self, cube, coord): """ - Returns a CF-netCDF variable name for the given coordinate. + 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. - * coord (:class:`iris.coords.Coord`): - An instance of a coordinate for which a CF-netCDF variable - name is required. + * coord (:class:`iris.coords.DimensionalMetadata`): + An instance of a coordinate (or similar), for which a CF-netCDF + variable name is required. Returns: A CF-netCDF variable name as a string. @@ -1812,6 +1810,8 @@ def _get_coord_variable_name(self, cube, coord): else: name = coord.standard_name or coord.long_name if not name or set(name).intersection(string.whitespace): + # NB don't know how to fix this if it is not a Coord or similar + assert isinstance(coord, iris.coords.Coord) # Auto-generate name based on associated dimensions. name = "" for dim in cube.coord_dims(coord): @@ -1825,166 +1825,57 @@ def _get_coord_variable_name(self, cube, coord): cf_name = self.cf_valid_var_name(cf_name) return cf_name - def _inner_create_cf_cellmeasure_or_ancil_variable( - self, cube, dimension_names, dimensional_metadata - ): - """ - Create the associated CF-netCDF variable in the netCDF dataset for the - given dimensional_metadata. - - Args: - - * cube (:class:`iris.cube.Cube`): - The associated cube being saved to CF-netCDF file. - * dimension_names (list): - Names for each dimension of the cube. - * dimensional_metadata (:class:`iris.coords.CellMeasure`): - A cell measure OR ancillary variable to be saved to the - CF-netCDF file. - In either case, provides data, units and standard/long/var names. - - Returns: - The string name of the associated CF-netCDF variable saved. - - """ - cf_name = self._get_coord_variable_name(cube, dimensional_metadata) - while cf_name in self._dataset.variables: - cf_name = self._increment_name(cf_name) - - # Derive the data dimension names for the coordinate. - cf_dimensions = [ - dimension_names[dim] - for dim in dimensional_metadata.cube_dims(cube) - ] - - # Get the data values. - data = dimensional_metadata.data - - if isinstance(dimensional_metadata, iris.coords.CellMeasure): - # Disallow saving of *masked* cell measures. - # NOTE: currently, this is the only functional difference required - # between variable creation for an ancillary and a cell measure. - if ma.is_masked(data): - # We can't save masked points properly, as we don't maintain a - # suitable fill_value. (Load will not record one, either). - msg = "Cell measures with missing data are not supported." - raise ValueError(msg) - - # Get the values in a form which is valid for the file format. - data = self._ensure_valid_dtype( - data, "coordinate", dimensional_metadata - ) - - # Create the CF-netCDF variable. - cf_var = self._dataset.createVariable( - cf_name, data.dtype.newbyteorder("="), cf_dimensions - ) - - # Add the data to the CF-netCDF variable. - cf_var[:] = data - - if dimensional_metadata.units.is_udunits(): - _setncattr(cf_var, "units", str(dimensional_metadata.units)) - - if dimensional_metadata.standard_name is not None: - _setncattr( - cf_var, "standard_name", dimensional_metadata.standard_name - ) - - if dimensional_metadata.long_name is not None: - _setncattr(cf_var, "long_name", dimensional_metadata.long_name) - - # Add any other custom coordinate attributes. - for name in sorted(dimensional_metadata.attributes): - value = dimensional_metadata.attributes[name] - - # Don't clobber existing attributes. - if not hasattr(cf_var, name): - _setncattr(cf_var, name, value) - - return cf_name - - def _create_cf_cell_measure_variable( - self, cube, dimension_names, cell_measure - ): - """ - Create the associated CF-netCDF variable in the netCDF dataset for the - given cell_measure. - - Args: - - * cube (:class:`iris.cube.Cube`): - The associated cube being saved to CF-netCDF file. - * dimension_names (list): - Names for each dimension of the cube. - * cell_measure (:class:`iris.coords.CellMeasure`): - The cell measure to be saved to CF-netCDF file. - - Returns: - The string name of the associated CF-netCDF variable saved. - - """ - # Note: currently shares variable creation code with ancillary-variables. - return self._inner_create_cf_cellmeasure_or_ancil_variable( - cube, dimension_names, cell_measure - ) + _ELEMENT_TYPE_NAMES = { + iris.coords.DimCoord: "coordinate", + iris.coords.AuxCoord: "coordinate", + iris.coords.CellMeasure: "cell-measure", + iris.coords.AncillaryVariable: "ancillary-variable", + } - def _create_cf_ancildata_variable( - self, cube, dimension_names, ancillary_variable - ): + def _create_generic_cf_array_var(self, cube, cube_dim_names, element): """ Create the associated CF-netCDF variable in the netCDF dataset for the - given ancillary variable. - - Args: - - * cube (:class:`iris.cube.Cube`): - The associated cube being saved to CF-netCDF file. - * dimension_names (list): - Names for each dimension of the cube. - * ancillary_variable (:class:`iris.coords.AncillaryVariable`): - The ancillary variable to be saved to the CF-netCDF file. - - Returns: - The string name of the associated CF-netCDF variable saved. + given dimensional_metadata. - """ - # Note: currently shares variable creation code with cell-measures. - return self._inner_create_cf_cellmeasure_or_ancil_variable( - cube, dimension_names, ancillary_variable - ) - - def _create_cf_coord_variable(self, cube, dimension_names, coord): - """ - Create the associated CF-netCDF variable in the netCDF dataset for the - given coordinate. If required, also create the CF-netCDF bounds - variable and associated dimension. + ..note:: + If the metadata element is a coord, it may also contain bounds. + In which case, an additional var is created and linked to it. Args: * cube (:class:`iris.cube.Cube`): The associated cube being saved to CF-netCDF file. - * dimension_names (list): - Names for each dimension of the cube. - * coord (:class:`iris.coords.Coord`): - The coordinate to be saved to CF-netCDF file. + * cube_dim_names (list of string): + The name of each dimension of the cube. + * element: + An Iris :class:`iris.coords.DimensionalMetadata`, belonging to the + cube. Provides data, units and standard/long/var names. Returns: - The string name of the associated CF-netCDF variable saved. + var_name (string): + The name of the CF-netCDF variable created. """ - cf_name = self._get_coord_variable_name(cube, coord) + # Work out the var-name to use. + cf_name = self._get_coord_variable_name(cube, element) while cf_name in self._dataset.variables: cf_name = self._increment_name(cf_name) - # Derive the data dimension names for the coordinate. - cf_dimensions = [ - dimension_names[dim] for dim in cube.coord_dims(coord) - ] - - if np.issubdtype(coord.points.dtype, np.str_): - string_dimension_depth = coord.points.dtype.itemsize - if coord.points.dtype.kind == "U": + # Get the list of file-dimensions (names), to create the variable. + element_dims = [ + cube_dim_names[dim] for dim in element.cube_dims(cube) + ] # NB using 'cube_dims' as this works for any type of element + + # Get the data values, in a way which works for any element type, as + # all are subclasses of _DimensionalMetadata. + # (e.g. =points if a coord, =data if an ancillary, etc) + data = element._values + + if np.issubdtype(data.dtype, np.str_): + # Deal with string-type variables. + # Typically CF label variables, but also possibly ancil-vars ? + string_dimension_depth = data.dtype.itemsize + if data.dtype.kind == "U": string_dimension_depth //= 4 string_dimension_name = "string%d" % string_dimension_depth @@ -1994,77 +1885,97 @@ def _create_cf_coord_variable(self, cube, dimension_names, coord): string_dimension_name, string_dimension_depth ) - # Add the string length dimension to dimension names. - cf_dimensions.append(string_dimension_name) + # Add the string length dimension to the variable dimensions. + element_dims.append(string_dimension_name) # Create the label coordinate variable. - cf_var = self._dataset.createVariable( - cf_name, "|S1", cf_dimensions - ) + cf_var = self._dataset.createVariable(cf_name, "|S1", element_dims) - # Add the payload to the label coordinate variable. - if len(cf_dimensions) == 1: - cf_var[:] = list( - "%- *s" % (string_dimension_depth, coord.points[0]) - ) + # Convert data from an array of strings into a character array + # with an extra string-length dimension. + if len(element_dims) == 1: + data = list("%- *s" % (string_dimension_depth, data[0])) else: - for index in np.ndindex(coord.points.shape): + orig_shape = data.shape + new_shape = orig_shape + (string_dimension_depth,) + new_data = np.zeros(new_shape, cf_var.dtype) + for index in np.ndindex(orig_shape): index_slice = tuple(list(index) + [slice(None, None)]) - cf_var[index_slice] = list( - "%- *s" % (string_dimension_depth, coord.points[index]) + new_data[index_slice] = list( + "%- *s" % (string_dimension_depth, data[index]) ) + data = new_data else: - # Identify the collection of coordinates that represent CF-netCDF - # coordinate variables. - cf_coordinates = cube.dim_coords + # A normal (numeric) variable. + # ensure a valid datatype for the file format. + element_type = self._ELEMENT_TYPE_NAMES[type(element)] + data = self._ensure_valid_dtype(data, element_type, element) + + # Check if this is a dim-coord. + is_dimcoord = element in cube.dim_coords + + if isinstance(element, iris.coords.CellMeasure): + # Disallow saving of *masked* cell measures. + # NOTE: currently, this is the only functional difference in + # variable creation between an ancillary and a cell measure. + if ma.is_masked(data): + # We can't save masked points properly, as we don't maintain + # a fill_value. (Load will not record one, either). + msg = "Cell measures with missing data are not supported." + raise ValueError(msg) - if coord in cf_coordinates: + if is_dimcoord: # By definition of a CF-netCDF coordinate variable this # coordinate must be 1-D and the name of the CF-netCDF variable # must be the same as its dimension name. - cf_name = cf_dimensions[0] - - # Get the values in a form which is valid for the file format. - points = self._ensure_valid_dtype( - coord.points, "coordinate", coord - ) + cf_name = element_dims[0] # Create the CF-netCDF variable. cf_var = self._dataset.createVariable( - cf_name, points.dtype.newbyteorder("="), cf_dimensions + cf_name, data.dtype.newbyteorder("="), element_dims ) # Add the axis attribute for spatio-temporal CF-netCDF coordinates. - if coord in cf_coordinates: - axis = iris.util.guess_coord_axis(coord) + if is_dimcoord: + axis = iris.util.guess_coord_axis(element) if axis is not None and axis.lower() in SPATIO_TEMPORAL_AXES: _setncattr(cf_var, "axis", axis.upper()) # Add the data to the CF-netCDF variable. - cf_var[:] = points + cf_var[:] = data # TODO: support dask streaming - # Create the associated CF-netCDF bounds variable. - self._create_cf_bounds(coord, cf_var, cf_name) + # Create the associated CF-netCDF bounds variable, if any. + self._create_cf_bounds(element, cf_var, cf_name) + + # Add the data to the CF-netCDF variable. + cf_var[:] = data # Deal with CF-netCDF units and standard name. - standard_name, long_name, units = self._cf_coord_identity(coord) + if isinstance(element, iris.coords.Coord): + # Fix "degree" units if needed. + # TODO: rewrite the handler routine to a more sensible API. + _, _, units_str = self._cf_coord_identity(element) + else: + units_str = str(element.units) - if cf_units.as_unit(units).is_udunits(): - _setncattr(cf_var, "units", units) + if cf_units.as_unit(units_str).is_udunits(): + _setncattr(cf_var, "units", units_str) + standard_name = element.standard_name if standard_name is not None: _setncattr(cf_var, "standard_name", standard_name) + long_name = element.long_name if long_name is not None: _setncattr(cf_var, "long_name", long_name) # Add the CF-netCDF calendar attribute. - if coord.units.calendar: - _setncattr(cf_var, "calendar", coord.units.calendar) + if element.units.calendar: + _setncattr(cf_var, "calendar", str(element.units.calendar)) # Add any other custom coordinate attributes. - for name in sorted(coord.attributes): - value = coord.attributes[name] + for name in sorted(element.attributes): + value = element.attributes[name] if name == "STASH": # Adopting provisional Metadata Conventions for representing MO diff --git a/lib/iris/tests/unit/fileformats/netcdf/test_Saver.py b/lib/iris/tests/unit/fileformats/netcdf/test_Saver.py index 2b0372dfa9..447cdd83eb 100644 --- a/lib/iris/tests/unit/fileformats/netcdf/test_Saver.py +++ b/lib/iris/tests/unit/fileformats/netcdf/test_Saver.py @@ -1031,7 +1031,7 @@ def test_masked_data__insitu(self): with self.temp_filename(".nc") as nc_path: saver = Saver(nc_path, "NETCDF4") with self.assertRaisesRegex(ValueError, self.exp_emsg): - saver._create_cf_cell_measure_variable( + saver._create_generic_cf_array_var( self.cube, self.names_map, self.cm ) From 69fd94ebdcdff39e5985b49d685f9e3b145bdd60 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 14 Sep 2021 16:45:08 +0100 Subject: [PATCH 02/26] Review changes. --- lib/iris/fileformats/netcdf.py | 79 ++++++++++++++++++---------------- 1 file changed, 42 insertions(+), 37 deletions(-) diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py index af2379a339..5be8f7c5ff 100644 --- a/lib/iris/fileformats/netcdf.py +++ b/lib/iris/fileformats/netcdf.py @@ -40,6 +40,7 @@ import iris.config import iris.coord_systems import iris.coords +from iris.coords import AncillaryVariable, AuxCoord, CellMeasure, DimCoord import iris.exceptions import iris.fileformats.cf import iris.io @@ -1356,40 +1357,50 @@ def _create_cf_dimensions( self._dataset.createDimension(dim_name, size) def _add_inner_related_vars( - self, - cube, - cf_var_cube, - dimension_names, - coordlike_elements, - role_attribute_name, + self, cube, cf_var_cube, dimension_names, coordlike_elements ): - # Common method to create a set of file variables and attach them to - # the parent data variable. - element_names = [] - # Add CF-netCDF variables for the associated auxiliary coordinates. - for element in sorted( - coordlike_elements, key=lambda element: element.name() - ): - # Create the associated CF-netCDF variable. - if element not in self._name_coord_map.coords: - cf_name = self._create_generic_cf_array_var( - cube, dimension_names, element - ) - self._name_coord_map.append(cf_name, element) + """ + Create a set of variables for aux-coords, ancillaries or cell-measures, + and attach them to the parent data variable. + + """ + if coordlike_elements: + # Choose the approriate parent attribute + elem_type = type(coordlike_elements[0]) + if elem_type in (AuxCoord, DimCoord): + role_attribute_name = "coordinates" + elif elem_type == AncillaryVariable: + role_attribute_name = "ancillary_variables" else: - cf_name = self._name_coord_map.name(element) + # We *only* handle aux-coords, cell-measures and ancillaries + assert elem_type == CellMeasure + role_attribute_name = "cell_measures" + + # Add CF-netCDF variables for the given cube components. + element_names = [] + for element in sorted( + coordlike_elements, key=lambda element: element.name() + ): + # Create the associated CF-netCDF variable. + if element not in self._name_coord_map.coords: + cf_name = self._create_generic_cf_array_var( + cube, dimension_names, element + ) + self._name_coord_map.append(cf_name, element) + else: + cf_name = self._name_coord_map.name(element) - if cf_name is not None: - if role_attribute_name == "cell_measures": - # In the case of cell-measures, the attribute entries are not just - # a var_name, but each have the form ": ". - cf_name = "{}: {}".format(element.measure, cf_name) - element_names.append(cf_name) + if cf_name is not None: + if role_attribute_name == "cell_measures": + # In the case of cell-measures, the attribute entries are not just + # a var_name, but each have the form ": ". + cf_name = "{}: {}".format(element.measure, cf_name) + element_names.append(cf_name) - # Add CF-netCDF references to the primary data variable. - if element_names: - variable_names = " ".join(sorted(element_names)) - _setncattr(cf_var_cube, role_attribute_name, variable_names) + # Add CF-netCDF references to the primary data variable. + if element_names: + variable_names = " ".join(sorted(element_names)) + _setncattr(cf_var_cube, role_attribute_name, variable_names) def _add_aux_coords(self, cube, cf_var_cube, dimension_names): """ @@ -1410,7 +1421,6 @@ def _add_aux_coords(self, cube, cf_var_cube, dimension_names): cf_var_cube, dimension_names, cube.aux_coords, - "coordinates", ) def _add_cell_measures(self, cube, cf_var_cube, dimension_names): @@ -1432,7 +1442,6 @@ def _add_cell_measures(self, cube, cf_var_cube, dimension_names): cf_var_cube, dimension_names, cube.cell_measures(), - "cell_measures", ) def _add_ancillary_variables(self, cube, cf_var_cube, dimension_names): @@ -1455,7 +1464,6 @@ def _add_ancillary_variables(self, cube, cf_var_cube, dimension_names): cf_var_cube, dimension_names, cube.ancillary_variables(), - "ancillary_variables", ) def _add_dim_coords(self, cube, dimension_names): @@ -1941,14 +1949,11 @@ def _create_generic_cf_array_var(self, cube, cube_dim_names, element): if axis is not None and axis.lower() in SPATIO_TEMPORAL_AXES: _setncattr(cf_var, "axis", axis.upper()) - # Add the data to the CF-netCDF variable. - cf_var[:] = data # TODO: support dask streaming - # Create the associated CF-netCDF bounds variable, if any. self._create_cf_bounds(element, cf_var, cf_name) # Add the data to the CF-netCDF variable. - cf_var[:] = data + cf_var[:] = data # TODO: support dask streaming # Deal with CF-netCDF units and standard name. if isinstance(element, iris.coords.Coord): From 306f6df507d465ccf329c3168a5b47bb4687c83a Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Wed, 25 Aug 2021 14:21:56 +0100 Subject: [PATCH 03/26] Ugrid save - first working. --- lib/iris/fileformats/netcdf.py | 346 +++++++++++++++--- .../experimental/test_ugrid_save.py | 96 +++++ 2 files changed, 382 insertions(+), 60 deletions(-) create mode 100644 lib/iris/tests/integration/experimental/test_ugrid_save.py diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py index 5be8f7c5ff..f79dda5ee8 100644 --- a/lib/iris/fileformats/netcdf.py +++ b/lib/iris/fileformats/netcdf.py @@ -977,6 +977,11 @@ def __setitem__(self, keys, arr): self.target[keys] = arr +# NOTE : this matches :class:`iris.experimental.ugrid.Mesh.LOCATIONS`, +# but in a specific preferred order. +MESH_LOCATIONS = ("face", "edge", "node") + + class Saver: """A manager for saving netcdf files.""" @@ -1192,15 +1197,20 @@ def write( self.check_attribute_compliance(coord, coord.points) # Get suitable dimension names. - dimension_names = self._get_dim_names(cube) + cube_dimensions, mesh_dimensions = self._get_dim_names(cube) - # Create the CF-netCDF data dimensions. - self._create_cf_dimensions(cube, dimension_names, unlimited_dimensions) + # Create all the CF-netCDF data dimensions. + # Put mesh dims first, then non-mesh dims in cube-occurring order. + nonmesh_dimensions = [ + dim for dim in cube_dimensions if dim not in mesh_dimensions + ] + all_dimensions = mesh_dimensions + nonmesh_dimensions + self._create_cf_dimensions(cube, all_dimensions, unlimited_dimensions) # Create the associated cube CF-netCDF data variable. cf_var_cube = self._create_cf_data_variable( cube, - dimension_names, + cube_dimensions, local_keys, zlib=zlib, complevel=complevel, @@ -1214,24 +1224,27 @@ def write( fill_value=fill_value, ) + # Add mesh components. + self._add_mesh(cube, cf_var_cube, cube_dimensions) + # Add coordinate variables. - self._add_dim_coords(cube, dimension_names) + self._add_dim_coords(cube, cube_dimensions) # Add the auxiliary coordinate variables and associate the data # variable to them - self._add_aux_coords(cube, cf_var_cube, dimension_names) + self._add_aux_coords(cube, cf_var_cube, cube_dimensions) # Add the cell_measures variables and associate the data # variable to them - self._add_cell_measures(cube, cf_var_cube, dimension_names) + self._add_cell_measures(cube, cf_var_cube, cube_dimensions) # Add the ancillary_variables variables and associate the data variable # to them - self._add_ancillary_variables(cube, cf_var_cube, dimension_names) + self._add_ancillary_variables(cube, cf_var_cube, cube_dimensions) # Add the formula terms to the appropriate cf variables for each # aux factory in the cube. - self._add_aux_factories(cube, cf_var_cube, dimension_names) + self._add_aux_factories(cube, cf_var_cube, cube_dimensions) # Add data variable-only attribute names to local_keys. if local_keys is None: @@ -1356,6 +1369,33 @@ def _create_cf_dimensions( size = self._existing_dim[dim_name] self._dataset.createDimension(dim_name, size) + def _add_mesh(self, cube, cf_var_cube, dimension_names): + """ + Add the cube's mesh, and all related variables to the dataset, + and associate it with the data variable. + Includes all the mesh-element coordinate and connectivity variables. + + Args: + + * cube (:class:`iris.cube.Cube`): + A :class:`iris.cube.Cube` to be saved to a netCDF file. + * cf_var_cube (:class:`netcdf.netcdf_variable`): + cf variable cube representation. + * dimension_names (list): + Names associated with the dimensions of the cube. + + """ + mesh = cube.mesh + if mesh: + if mesh in self._name_coord_map.coords: + cf_name = self._name_coord_map.name(mesh) + else: + cf_name = self._create_mesh(cube, dimension_names, mesh) + self._name_coord_map.append(cf_name, mesh) + + _setncattr(cf_var_cube, "mesh", cf_name) + _setncattr(cf_var_cube, "location", cube.location) + def _add_inner_related_vars( self, cube, cf_var_cube, dimension_names, coordlike_elements ): @@ -1416,11 +1456,15 @@ def _add_aux_coords(self, cube, cf_var_cube, dimension_names): Names associated with the dimensions of the cube. """ + # Exclude mesh coords, which are bundled in with the aux-coords. + aux_coords_no_mesh = [ + coord for coord in cube.aux_coords if not hasattr(coord, "mesh") + ] return self._add_inner_related_vars( cube, cf_var_cube, dimension_names, - cube.aux_coords, + aux_coords_no_mesh, ) def _add_cell_measures(self, cube, cf_var_cube, dimension_names): @@ -1593,39 +1637,94 @@ def _get_dim_names(self, cube): A :class:`iris.cube.Cube` to be saved to a netCDF file. Returns: - List of dimension names with length equal the number of dimensions - in the cube. + cube_dimensions, mesh_dimensions + * cube_dimensions (list of string): + A lists of dimension names for each dimension of the cube + * mesh_dimensions (list of string): + A list of the mesh dimensions of the attached mesh, if any. + + ..note:: + One of the mesh dimensions will generally also appear in the cube + dimensions. """ dimension_names = [] + + def record_dimension(dim_name, length, dim_coords=[]): + """ + Record a file dimension and its length. + + If the dimension has been seen already, check that it's length + matches the earlier finding. + + """ + if dim_name not in self._existing_dim: + self._existing_dim[dim_name] = length + else: + # Just make sure we never re-write one. + # TODO: possibly merits a proper Exception with message + assert self._existing_dim[dim_name] == length + + # Add new coords (mesh or dim) to the already-seen list. + for coord in dim_coords: + if coord not in self._dim_coords: + self._dim_coords.append(coord) + + # Add the latest name to the returned list + dimension_names.append(dim_name) + + # Get info on mesh, first. + dimension_names.clear() + mesh = cube.mesh + if mesh is not None: + # Create all the mesh dimensions. + # NOTE: one of these will be a cube dimension, but that is not + # special. We *do* want to list these in a priority order + # (face,edge,node), and before non-mesh dimensions. + mesh_coords = cube.coords(mesh_coords=True) + + for location in MESH_LOCATIONS: + location_select_key = f"include_{location}s" + mesh_coords = mesh.coords(**{location_select_key: True}) + if len(mesh_coords) > 0: + (dim_length,) = mesh_coords[0].shape # Must be 1-d + location_dim_attr = f"{location}_dimension" + dim_name = getattr(mesh, location_dim_attr) + if dim_name is None: + dim_name = f"{mesh.name()}_{location}" + while dim_name in self._existing_dim: + self._increment_name(dim_name) + + record_dimension(dim_name, dim_length) + + mesh_dimensions = dimension_names.copy() + + # Get the cube dimensions, in order. + dimension_names.clear() for dim in range(cube.ndim): - coords = cube.coords(dimensions=dim, dim_coords=True) - if coords: - coord = coords[0] + # 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 - dim_name = self._get_coord_variable_name(cube, coord) # Add only dimensions that have not already been added. - if coord not in self._dim_coords: - # Determine unique dimension name + if coord in self._dim_coords: + # Return the dim_name associated with the existing + # coordinate. + dim_name = self._name_coord_map.name(coord) + else: + # Determine 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) - # Update names added, current cube dim names used and - # unique coordinates added. - self._existing_dim[dim_name] = coord.shape[0] - dimension_names.append(dim_name) - self._dim_coords.append(coord) - else: - # Return the dim_name associated with the existing - # coordinate. - dim_name = self._name_coord_map.name(coord) - dimension_names.append(dim_name) - else: # No CF-netCDF coordinates describe this data dimension. + # Make up a new, distinct dimension name dim_name = "dim%d" % dim if dim_name in self._existing_dim: # Increment name if conflicted with one already existing. @@ -1636,14 +1735,13 @@ def _get_dim_names(self, cube): or dim_name in self._name_coord_map.names ): dim_name = self._increment_name(dim_name) - # Update dictionary with new entry - self._existing_dim[dim_name] = cube.shape[dim] - else: - # Update dictionary with new entry - self._existing_dim[dim_name] = cube.shape[dim] - dimension_names.append(dim_name) - return dimension_names + # Record the dimension. + record_dimension(dim_name, cube.shape[dim], dim_coords) + + cube_dimensions = dimension_names.copy() + + return cube_dimensions, mesh_dimensions @staticmethod def cf_valid_var_name(var_name): @@ -1818,29 +1916,149 @@ def _get_coord_variable_name(self, cube, coord): else: name = coord.standard_name or coord.long_name if not name or set(name).intersection(string.whitespace): - # NB don't know how to fix this if it is not a Coord or similar - assert isinstance(coord, iris.coords.Coord) - # Auto-generate name based on associated dimensions. - name = "" - for dim in cube.coord_dims(coord): - name += "dim{}".format(dim) - # Handle scalar coordinate (dims == ()). - if not name: - name = "unknown_scalar" + # We need to invent a name, based on its associated dimensions. + if cube.coords(coord): + # It is a regular cube coordinate. + # Auto-generate a name based on the dims. + name = "" + for dim in cube.coord_dims(coord): + name += "dim{}".format(dim) + # Handle scalar coordinate (dims == ()). + if not name: + name = "unknown_scalar" + else: + # Not a cube coord, so must be a connectivity or + # element-coordinate of the mesh. + # Name it for it's first dim, i.e. mesh-dim of its location. + from iris.experimental.ugrid import Connectivity + + if isinstance(coord, Connectivity): + location = coord.cf_role.split("_")[0] + else: + # Must be a mesh-element coordinate. + location = None + for test_location in MESH_LOCATIONS: + include_key = f"include_{test_location}s" + if coord in cube.mesh.coords({include_key: True}): + location = test_location + break + assert location is not None + location_dim_attr = f"{location}_dimension" + name = getattr(cube.mesh, location_dim_attr) + # Convert to lower case and replace whitespace by underscores. cf_name = "_".join(name.lower().split()) cf_name = self.cf_valid_var_name(cf_name) return cf_name - _ELEMENT_TYPE_NAMES = { - iris.coords.DimCoord: "coordinate", - iris.coords.AuxCoord: "coordinate", - iris.coords.CellMeasure: "cell-measure", - iris.coords.AncillaryVariable: "ancillary-variable", - } + def _get_mesh_variable_name(self, mesh): + """ + Returns a CF-netCDF variable name for the given coordinate. + + Args: + + * mesh (:class:`iris.experimental.ugrid.Mesh`): + An instance of a Mesh for which a CF-netCDF variable name is + required. + + Returns: + A CF-netCDF variable name as a string. + + """ + cf_name = mesh.var_name + if mesh.var_name is None: + # Prefer var-name, but accept long/standard as aliases. + cf_name = mesh.var_name or mesh.long_name or mesh.standard_name + if not cf_name: + # Auto-generate a name based on mesh properties. + cf_name = f"Mesh_{mesh.topology_dimension}d" + + # Ensure valid form for var-name. + cf_name = self.cf_valid_var_name(cf_name) + return cf_name - def _create_generic_cf_array_var(self, cube, cube_dim_names, element): + def _create_mesh(self, cube, dimension_names, mesh): + """ + Create all the variables associated with a mesh in the netCDF dataset. + + Args: + + * cube (:class:`iris.cube.Cube`): + The associated cube being saved to CF-netCDF file. + * dimension_names (list): + Names for each dimension of the cube. + * mesh (:class:`iris.experimental.ugrid.Mesh`): + The Mesh to be saved to CF-netCDF file. + + Returns: + The string name of the associated CF-netCDF variable saved. + + """ + # First choose a var-name for the mesh variable itself. + cf_mesh_name = self._get_mesh_variable_name(mesh) + # Disambiguate any possible clashes. + while cf_mesh_name in self._dataset.variables: + cf_mesh_name = self._increment_name(cf_mesh_name) + + # Create the main variable + cf_mesh_var = self._dataset.createVariable( + cf_mesh_name, + np.dtype(np.int64), + [], + ) + + # Get the mesh-element dim names. + mesh_dims = { + location: getattr(mesh, f"{location}_dimension") + for location in MESH_LOCATIONS + } + # Add the element coordinate variables. + for location in MESH_LOCATIONS: + coords_attr_name = f"{location}_coords" + mesh_coords = getattr(mesh, coords_attr_name) + if mesh_coords: + coord_names = [] + for coord in mesh_coords: + if coord is None: + continue # an awkward thing that mesh.coords does + coord_name = self._create_generic_cf_array_var( + cube, [], coord, element_dims=(mesh_dims[location],) + ) + coord_names.append(coord_name) + # Record the coordinates (if any) on the mesh variable. + if coord_names: + coord_names = " ".join(coord_names) + _setncattr(cf_mesh_var, coords_attr_name, coord_names) + + # Add all the connectivity variables. + # pre-fetch the set + ignore "None"s -- looks like a bug ? + conns = [conn for conn in mesh.all_connectivities if conn is not None] + for conn in conns: + # Get the connectivity role, = "{loc1}_{loc2}_connectivity". + cf_conn_attr_name = conn.cf_role + loc_from, loc_to, _ = cf_conn_attr_name.split("_") + # Construct a trailing dimension name. + last_dim = f"{cf_mesh_name}_N_{loc_from}_{loc_to}s" + # Create if it does not already exist. + if last_dim not in self._dataset.dimensions: + length = conn.shape[1] + self._dataset.createDimension(last_dim, length) + # Create variable. + conn_dims = (mesh_dims[loc_from], last_dim) + cf_conn_name = self._create_generic_cf_array_var( + cube, [], conn, element_dims=conn_dims + ) + # Record in a mesh var attr whose name is the cf_role. + _setncattr(cf_mesh_var, cf_conn_attr_name, cf_conn_name) + + # TODO: possibly we need to handle boundary coords/conns separately + + return cf_mesh_name + + def _create_generic_cf_array_var( + self, cube, cube_dim_names, element, element_dims=None + ): """ Create the associated CF-netCDF variable in the netCDF dataset for the given dimensional_metadata. @@ -1855,9 +2073,16 @@ def _create_generic_cf_array_var(self, cube, cube_dim_names, element): The associated cube being saved to CF-netCDF file. * cube_dim_names (list of string): The name of each dimension of the cube. - * element: - An Iris :class:`iris.coords.DimensionalMetadata`, belonging to the - cube. Provides data, units and standard/long/var names. + Not used if 'element_dims' is not None. + * element (:class:`iris.coords.DimensionalMetadata`): + A cube component, represented by a file variable, e.g. a coordinate + or cell-measure. + Provides data, units and standard/long/var names. + * element_dims (list of string, or None): + If set, contains the variable dimension (names), + otherwise these are taken from `element.cube_dims[cube]`. + For Mesh components (element coordinates and connectivities), this + *must* be passed in, as "element.cube_dims" does not function. Returns: var_name (string): @@ -1869,10 +2094,11 @@ def _create_generic_cf_array_var(self, cube, cube_dim_names, element): while cf_name in self._dataset.variables: cf_name = self._increment_name(cf_name) - # Get the list of file-dimensions (names), to create the variable. - element_dims = [ - cube_dim_names[dim] for dim in element.cube_dims(cube) - ] # NB using 'cube_dims' as this works for any type of element + if element_dims is None: + # Get the list of file-dimensions (names), to create the variable. + element_dims = [ + cube_dim_names[dim] for dim in element.cube_dims(cube) + ] # NB using 'cube_dims' as this works for any type of element # Get the data values, in a way which works for any element type, as # all are subclasses of _DimensionalMetadata. @@ -1916,7 +2142,7 @@ def _create_generic_cf_array_var(self, cube, cube_dim_names, element): else: # A normal (numeric) variable. # ensure a valid datatype for the file format. - element_type = self._ELEMENT_TYPE_NAMES[type(element)] + element_type = type(element).__name__ data = self._ensure_valid_dtype(data, element_type, element) # Check if this is a dim-coord. diff --git a/lib/iris/tests/integration/experimental/test_ugrid_save.py b/lib/iris/tests/integration/experimental/test_ugrid_save.py new file mode 100644 index 0000000000..517e8d1cb3 --- /dev/null +++ b/lib/iris/tests/integration/experimental/test_ugrid_save.py @@ -0,0 +1,96 @@ +# Copyright Iris contributors +# +# This file is part of Iris and is released under the LGPL license. +# See COPYING and COPYING.LESSER in the root of the repository for full +# licensing details. +""" +Integration tests for NetCDF-UGRID file saving. + +""" +# Import iris.tests first so that some things can be initialised before +# importing anything else. +import iris.tests as tests # isort:skip + +from pathlib import Path +import shutil +import tempfile + +import iris +from iris.experimental.ugrid import PARSE_UGRID_ON_LOAD +import iris.fileformats.netcdf +from iris.tests import IrisTest +from iris.tests.stock.mesh import ( # sample_mesh,; sample_meshcoord, + sample_mesh_cube, +) + + +class TestBasicSave(IrisTest): + @classmethod + def setUpClass(cls): + cls.temp_dir = Path(tempfile.mkdtemp()) + + @classmethod + def tearDownClass(cls): + shutil.rmtree(cls.temp_dir) + + def test_parts(self): + cube = sample_mesh_cube() + # Make this REALLY minimal + # Note: removing the 'extras' in the 'sample+mesh_cube' + cube = cube[0] # Strip out first dimension. + # # Remove un-needed extra coords (for now). + for name in ("level", "i_mesh_face", "mesh_face_aux"): + cube.remove_coord(name) + print(cube) + + temp_path = str(self.temp_dir / "tmp.nc") + with iris.fileformats.netcdf.Saver(temp_path, "NETCDF4") as saver: + saver.write(cube) + # TODO: do some actual testing, beyond "does not fail" + # self.assertCDL(temp_path) # TODO: for now onl + + def test_basic_save(self): + # Generate an ultra-simple unstructured cube + cube = sample_mesh_cube() + print(cube) + return + + # Save it out, and check that the CDL is as expected. + tempfile_path = str(self.temp_dir / "basic.nc") + iris.save(cube, tempfile_path) + self.assertCDL(tempfile_path) + + # def test_complex_multiple_save(self): + # source_dirpath = Path(iris.tests.get_data_path( + # ['NetCDF', 'unstructured_grid', 'lfric_surface_mean.nc'])) + # + # # Save it out, and check that the CDL is as expected. + # tempfile_path = str(self.temp_dir / 'basic.nc') + # iris.save(cube, tempfile_path) + # self.assertCDL(tempfile_path) + + def test_roundtrip(self): + source_dirpath = Path( + iris.tests.get_data_path(["NetCDF", "unstructured_grid"]) + ) + file_name = "data_C4.nc" + source_filepath = str(source_dirpath / file_name) + + # # Ensure that our test reference text matches the original input file + # self.assertCDL(source_filepath) + + # Load the data (one cube) and save it out to a temporary file of the same name + with PARSE_UGRID_ON_LOAD.context(): + cube = iris.load_cube(source_filepath) + + print(cube) + target_filepath = str(self.temp_dir / file_name) + iris.save(cube, target_filepath) + + print(cube) + # # Ensure that the saved result is identical + # self.assertCDL(target_filepath) + + +if __name__ == "__main__": + tests.main() From cab7e34af7b55f7f004514b5a62ce031b84ff9b0 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 31 Aug 2021 19:20:49 +0100 Subject: [PATCH 04/26] Required mesh attributes, loadback check. --- lib/iris/fileformats/netcdf.py | 96 +++++++++++-------- .../experimental/test_ugrid_save.py | 29 +++++- 2 files changed, 84 insertions(+), 41 deletions(-) diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py index f79dda5ee8..277f4450b7 100644 --- a/lib/iris/fileformats/netcdf.py +++ b/lib/iris/fileformats/netcdf.py @@ -2008,6 +2008,12 @@ def _create_mesh(self, cube, dimension_names, mesh): [], ) + # Add the basic essential attributes + _setncattr(cf_mesh_var, "cf_role", "mesh_topology") + _setncattr(cf_mesh_var, "topology_dimension", mesh.topology_dimension) + # Add 'standard' names + units atributes + self._set_cf_var_attributes(cf_mesh_var, mesh) + # Get the mesh-element dim names. mesh_dims = { location: getattr(mesh, f"{location}_dimension") @@ -2015,8 +2021,9 @@ def _create_mesh(self, cube, dimension_names, mesh): } # Add the element coordinate variables. for location in MESH_LOCATIONS: - coords_attr_name = f"{location}_coords" - mesh_coords = getattr(mesh, coords_attr_name) + coords_meshobj_attr = f"{location}_coords" + coords_file_attr = f"{location}_coordinates" + mesh_coords = getattr(mesh, coords_meshobj_attr) if mesh_coords: coord_names = [] for coord in mesh_coords: @@ -2029,7 +2036,7 @@ def _create_mesh(self, cube, dimension_names, mesh): # Record the coordinates (if any) on the mesh variable. if coord_names: coord_names = " ".join(coord_names) - _setncattr(cf_mesh_var, coords_attr_name, coord_names) + _setncattr(cf_mesh_var, coords_file_attr, coord_names) # Add all the connectivity variables. # pre-fetch the set + ignore "None"s -- looks like a bug ? @@ -2049,13 +2056,56 @@ def _create_mesh(self, cube, dimension_names, mesh): cf_conn_name = self._create_generic_cf_array_var( cube, [], conn, element_dims=conn_dims ) - # Record in a mesh var attr whose name is the cf_role. + # Add essential attributes to the Connectivity variable. + cf_conn_var = self._dataset.variables[cf_conn_name] + _setncattr(cf_conn_var, "cf_role", cf_conn_attr_name) + _setncattr(cf_conn_var, "start_index", conn.start_index) + + # Record the connectivity on the parent mesh var. _setncattr(cf_mesh_var, cf_conn_attr_name, cf_conn_name) # TODO: possibly we need to handle boundary coords/conns separately return cf_mesh_name + def _set_cf_var_attributes(self, cf_var, element): + # Deal with CF-netCDF units and standard name. + if isinstance(element, iris.coords.Coord): + # Fix "degree" units if needed. + # TODO: rewrite the handler routine to a more sensible API. + _, _, units_str = self._cf_coord_identity(element) + else: + units_str = str(element.units) + + if cf_units.as_unit(units_str).is_udunits(): + _setncattr(cf_var, "units", units_str) + + standard_name = element.standard_name + if standard_name is not None: + _setncattr(cf_var, "standard_name", standard_name) + + long_name = element.long_name + if long_name is not None: + _setncattr(cf_var, "long_name", long_name) + + # Add the CF-netCDF calendar attribute. + if element.units.calendar: + _setncattr(cf_var, "calendar", str(element.units.calendar)) + + # Add any other custom coordinate attributes. + for name in sorted(element.attributes): + value = element.attributes[name] + + if name == "STASH": + # Adopting provisional Metadata Conventions for representing MO + # Scientific Data encoded in NetCDF Format. + name = "um_stash_source" + value = str(value) + + # Don't clobber existing attributes. + if not hasattr(cf_var, name): + _setncattr(cf_var, name, value) + def _create_generic_cf_array_var( self, cube, cube_dim_names, element, element_dims=None ): @@ -2181,42 +2231,8 @@ def _create_generic_cf_array_var( # Add the data to the CF-netCDF variable. cf_var[:] = data # TODO: support dask streaming - # Deal with CF-netCDF units and standard name. - if isinstance(element, iris.coords.Coord): - # Fix "degree" units if needed. - # TODO: rewrite the handler routine to a more sensible API. - _, _, units_str = self._cf_coord_identity(element) - else: - units_str = str(element.units) - - if cf_units.as_unit(units_str).is_udunits(): - _setncattr(cf_var, "units", units_str) - - standard_name = element.standard_name - if standard_name is not None: - _setncattr(cf_var, "standard_name", standard_name) - - long_name = element.long_name - if long_name is not None: - _setncattr(cf_var, "long_name", long_name) - - # Add the CF-netCDF calendar attribute. - if element.units.calendar: - _setncattr(cf_var, "calendar", str(element.units.calendar)) - - # Add any other custom coordinate attributes. - for name in sorted(element.attributes): - value = element.attributes[name] - - if name == "STASH": - # Adopting provisional Metadata Conventions for representing MO - # Scientific Data encoded in NetCDF Format. - name = "um_stash_source" - value = str(value) - - # Don't clobber existing attributes. - if not hasattr(cf_var, name): - _setncattr(cf_var, name, value) + # Add names + units + self._set_cf_var_attributes(cf_var, element) return cf_name diff --git a/lib/iris/tests/integration/experimental/test_ugrid_save.py b/lib/iris/tests/integration/experimental/test_ugrid_save.py index 517e8d1cb3..cdcb1beaf8 100644 --- a/lib/iris/tests/integration/experimental/test_ugrid_save.py +++ b/lib/iris/tests/integration/experimental/test_ugrid_save.py @@ -87,10 +87,37 @@ def test_roundtrip(self): target_filepath = str(self.temp_dir / file_name) iris.save(cube, target_filepath) - print(cube) + from subprocess import check_output + + print("") + print("ORIGINAL:") + text = check_output("ncdump -h " + source_filepath, shell=True) + print(text.decode()) + print("") + print("RE-SAVED:") + text = check_output("ncdump -h " + target_filepath, shell=True) + print(text.decode()) # # Ensure that the saved result is identical # self.assertCDL(target_filepath) + # Now try loading BACK... + with PARSE_UGRID_ON_LOAD.context(): + cube_reload = iris.load_cube(target_filepath) + print("") + print("OUTPUT-RE-LOADED:") + print(cube_reload) + + mesh_orig = cube.mesh + mesh_reload = cube_reload.mesh + for propname in dir(mesh_orig): + if not propname.startswith("_"): + prop_orig = getattr(mesh_orig, propname) + if not callable(prop_orig): + prop_reload = getattr(mesh_reload, propname) + self.assertEqual( + (propname, prop_reload), (propname, prop_orig) + ) + if __name__ == "__main__": tests.main() From 89864f058f0cd32e6c6ccf5838cc0b11d2dc09d6 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 2 Sep 2021 11:54:09 +0100 Subject: [PATCH 05/26] Odd fixes to netcdf-save. --- lib/iris/fileformats/netcdf.py | 133 ++++++++++++++++++++------------- 1 file changed, 81 insertions(+), 52 deletions(-) diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py index 277f4450b7..19f973191f 100644 --- a/lib/iris/fileformats/netcdf.py +++ b/lib/iris/fileformats/netcdf.py @@ -979,7 +979,7 @@ def __setitem__(self, keys, arr): # NOTE : this matches :class:`iris.experimental.ugrid.Mesh.LOCATIONS`, # but in a specific preferred order. -MESH_LOCATIONS = ("face", "edge", "node") +MESH_LOCATIONS = ("node", "edge", "face") class Saver: @@ -1207,6 +1207,11 @@ def write( all_dimensions = mesh_dimensions + nonmesh_dimensions self._create_cf_dimensions(cube, all_dimensions, unlimited_dimensions) + # Create the mesh components, if there is a mesh. + # We do this before creating the data-var, so that mesh vars precede + # data-vars in the file. + cf_mesh_name = self._add_mesh(cube) + # Create the associated cube CF-netCDF data variable. cf_var_cube = self._create_cf_data_variable( cube, @@ -1224,8 +1229,12 @@ def write( fill_value=fill_value, ) - # Add mesh components. - self._add_mesh(cube, cf_var_cube, cube_dimensions) + # Associate any mesh with the data-variable. + # N.B. _add_mesh cannot do this, as we want to put mesh variables + # before data-variables in the file. + if cf_mesh_name is not None: + _setncattr(cf_var_cube, "mesh", cf_mesh_name) + _setncattr(cf_var_cube, "location", cube.location) # Add coordinate variables. self._add_dim_coords(cube, cube_dimensions) @@ -1369,32 +1378,37 @@ def _create_cf_dimensions( size = self._existing_dim[dim_name] self._dataset.createDimension(dim_name, size) - def _add_mesh(self, cube, cf_var_cube, dimension_names): + def _add_mesh(self, cube): """ - Add the cube's mesh, and all related variables to the dataset, - and associate it with the data variable. + Add the cube's mesh, and all related variables to the dataset. Includes all the mesh-element coordinate and connectivity variables. + ..note:: + + Here, we do *not* add the relevant referencing attributes to the + data-variable, because we want to create the data-variable later. + Args: * cube (:class:`iris.cube.Cube`): A :class:`iris.cube.Cube` to be saved to a netCDF file. - * cf_var_cube (:class:`netcdf.netcdf_variable`): - cf variable cube representation. - * dimension_names (list): - Names associated with the dimensions of the cube. + + Returns: + * cf_mesh_name (string or None): + The name of the mesh variable created, or None if the cube does not + have a mesh. """ + cf_mesh_name = None mesh = cube.mesh if mesh: if mesh in self._name_coord_map.coords: - cf_name = self._name_coord_map.name(mesh) + cf_mesh_name = self._name_coord_map.name(mesh) else: - cf_name = self._create_mesh(cube, dimension_names, mesh) - self._name_coord_map.append(cf_name, mesh) + cf_mesh_name = self._create_mesh(cube, mesh) + self._name_coord_map.append(cf_mesh_name, mesh) - _setncattr(cf_var_cube, "mesh", cf_name) - _setncattr(cf_var_cube, "location", cube.location) + return cf_mesh_name def _add_inner_related_vars( self, cube, cf_var_cube, dimension_names, coordlike_elements @@ -1676,13 +1690,15 @@ def record_dimension(dim_name, length, dim_coords=[]): # Get info on mesh, first. dimension_names.clear() mesh = cube.mesh - if mesh is not None: + if mesh is None: + mesh_location_dimnames = {} + cube_mesh_dim = None + else: # Create all the mesh dimensions. # NOTE: one of these will be a cube dimension, but that is not # special. We *do* want to list these in a priority order # (face,edge,node), and before non-mesh dimensions. - mesh_coords = cube.coords(mesh_coords=True) - + mesh_location_dimnames = {} for location in MESH_LOCATIONS: location_select_key = f"include_{location}s" mesh_coords = mesh.coords(**{location_select_key: True}) @@ -1696,46 +1712,57 @@ def record_dimension(dim_name, length, dim_coords=[]): self._increment_name(dim_name) record_dimension(dim_name, dim_length) + mesh_location_dimnames[location] = dim_name + + # Identify the cube dimension which maps to the mesh. + # NB always at least one mesh-coord: it has exactly one cube-dim + (cube_mesh_dim,) = cube.coords(mesh_coords=True)[0].cube_dims(cube) mesh_dimensions = dimension_names.copy() # Get the cube dimensions, in order. dimension_names.clear() for dim in range(cube.ndim): - # 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. - if coord in self._dim_coords: - # Return the dim_name associated with the existing - # coordinate. - dim_name = self._name_coord_map.name(coord) - else: - # Determine 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 dim == cube_mesh_dim: + # Handle a mesh dimension: we already named this. + dim_coords = [] + dim_name = mesh_location_dimnames[cube.location] else: - # No CF-netCDF coordinates describe this data dimension. - # Make up a new, distinct dimension name - dim_name = "dim%d" % 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. + if coord in self._dim_coords: + # Return the dim_name associated with the existing + # coordinate. + dim_name = self._name_coord_map.name(coord) + else: + # Determine 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) + else: + # No CF-netCDF coordinates describe this data dimension. + # Make up a new, distinct dimension name + dim_name = "dim%d" % 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(dim_name, cube.shape[dim], dim_coords) @@ -1978,7 +2005,7 @@ def _get_mesh_variable_name(self, mesh): cf_name = self.cf_valid_var_name(cf_name) return cf_name - def _create_mesh(self, cube, dimension_names, mesh): + def _create_mesh(self, cube, mesh): """ Create all the variables associated with a mesh in the netCDF dataset. @@ -1986,8 +2013,6 @@ def _create_mesh(self, cube, dimension_names, mesh): * cube (:class:`iris.cube.Cube`): The associated cube being saved to CF-netCDF file. - * dimension_names (list): - Names for each dimension of the cube. * mesh (:class:`iris.experimental.ugrid.Mesh`): The Mesh to be saved to CF-netCDF file. @@ -2004,13 +2029,17 @@ def _create_mesh(self, cube, dimension_names, mesh): # Create the main variable cf_mesh_var = self._dataset.createVariable( cf_mesh_name, - np.dtype(np.int64), + np.dtype(np.int32), [], ) # Add the basic essential attributes _setncattr(cf_mesh_var, "cf_role", "mesh_topology") - _setncattr(cf_mesh_var, "topology_dimension", mesh.topology_dimension) + _setncattr( + cf_mesh_var, + "topology_dimension", + np.int32(mesh.topology_dimension), + ) # Add 'standard' names + units atributes self._set_cf_var_attributes(cf_mesh_var, mesh) @@ -2023,7 +2052,7 @@ def _create_mesh(self, cube, dimension_names, mesh): for location in MESH_LOCATIONS: coords_meshobj_attr = f"{location}_coords" coords_file_attr = f"{location}_coordinates" - mesh_coords = getattr(mesh, coords_meshobj_attr) + mesh_coords = getattr(mesh, coords_meshobj_attr, None) if mesh_coords: coord_names = [] for coord in mesh_coords: @@ -2046,7 +2075,7 @@ def _create_mesh(self, cube, dimension_names, mesh): cf_conn_attr_name = conn.cf_role loc_from, loc_to, _ = cf_conn_attr_name.split("_") # Construct a trailing dimension name. - last_dim = f"{cf_mesh_name}_N_{loc_from}_{loc_to}s" + last_dim = f"{cf_mesh_name}_{loc_from}_N_{loc_to}s" # Create if it does not already exist. if last_dim not in self._dataset.dimensions: length = conn.shape[1] From 0aab860c065e18f72f4683eed93fba28bc496719 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 2 Sep 2021 14:22:04 +0100 Subject: [PATCH 06/26] Added tests based on UGRID cdl examples. --- .../experimental/test_ugrid_save.py | 175 +++++++++--------- .../ugrid_conventions_examples/README.txt | 46 +++++ .../ugrid_ex1_1d_mesh.cdl | 55 ++++++ .../ugrid_ex2_2d_triangular.cdl | 84 +++++++++ .../ugrid_ex3_2d_flexible.cdl | 99 ++++++++++ .../ugrid_ex4_3d_layered.cdl | 124 +++++++++++++ .../TestBasicSave/ugrid_ex1_1d_mesh.cdl | 39 ++++ .../TestBasicSave/ugrid_ex2_2d_triangular.cdl | 75 ++++++++ .../TestBasicSave/ugrid_ex3_2d_flexible.cdl | 75 ++++++++ .../TestBasicSave/ugrid_ex4_3d_layered.cdl | 106 +++++++++++ lib/iris/tests/stock/netcdf.py | 8 +- 11 files changed, 797 insertions(+), 89 deletions(-) create mode 100644 lib/iris/tests/integration/experimental/ugrid_conventions_examples/README.txt create mode 100644 lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex1_1d_mesh.cdl create mode 100644 lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex2_2d_triangular.cdl create mode 100644 lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex3_2d_flexible.cdl create mode 100644 lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex4_3d_layered.cdl create mode 100644 lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex1_1d_mesh.cdl create mode 100644 lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex2_2d_triangular.cdl create mode 100644 lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex3_2d_flexible.cdl create mode 100644 lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex4_3d_layered.cdl diff --git a/lib/iris/tests/integration/experimental/test_ugrid_save.py b/lib/iris/tests/integration/experimental/test_ugrid_save.py index cdcb1beaf8..969b5d7557 100644 --- a/lib/iris/tests/integration/experimental/test_ugrid_save.py +++ b/lib/iris/tests/integration/experimental/test_ugrid_save.py @@ -11,112 +11,113 @@ # importing anything else. import iris.tests as tests # isort:skip +import glob from pathlib import Path import shutil +from subprocess import check_call import tempfile import iris from iris.experimental.ugrid import PARSE_UGRID_ON_LOAD import iris.fileformats.netcdf from iris.tests import IrisTest -from iris.tests.stock.mesh import ( # sample_mesh,; sample_meshcoord, - sample_mesh_cube, -) +from iris.tests.stock.netcdf import _add_standard_data class TestBasicSave(IrisTest): @classmethod def setUpClass(cls): cls.temp_dir = Path(tempfile.mkdtemp()) + cls.examples_dir = ( + Path(__file__).absolute().parent / "ugrid_conventions_examples" + ) + example_paths = glob.glob(str(cls.examples_dir / "*ex*.cdl")) + example_names = [ + str(Path(filepath).name).split("_")[1] # = "ex" + for filepath in example_paths + ] + cls.example_names_paths = { + name: path for name, path in zip(example_names, example_paths) + } @classmethod def tearDownClass(cls): shutil.rmtree(cls.temp_dir) - def test_parts(self): - cube = sample_mesh_cube() - # Make this REALLY minimal - # Note: removing the 'extras' in the 'sample+mesh_cube' - cube = cube[0] # Strip out first dimension. - # # Remove un-needed extra coords (for now). - for name in ("level", "i_mesh_face", "mesh_face_aux"): - cube.remove_coord(name) - print(cube) - - temp_path = str(self.temp_dir / "tmp.nc") - with iris.fileformats.netcdf.Saver(temp_path, "NETCDF4") as saver: - saver.write(cube) - # TODO: do some actual testing, beyond "does not fail" - # self.assertCDL(temp_path) # TODO: for now onl - - def test_basic_save(self): - # Generate an ultra-simple unstructured cube - cube = sample_mesh_cube() - print(cube) - return - - # Save it out, and check that the CDL is as expected. - tempfile_path = str(self.temp_dir / "basic.nc") - iris.save(cube, tempfile_path) - self.assertCDL(tempfile_path) - - # def test_complex_multiple_save(self): - # source_dirpath = Path(iris.tests.get_data_path( - # ['NetCDF', 'unstructured_grid', 'lfric_surface_mean.nc'])) - # - # # Save it out, and check that the CDL is as expected. - # tempfile_path = str(self.temp_dir / 'basic.nc') - # iris.save(cube, tempfile_path) - # self.assertCDL(tempfile_path) - - def test_roundtrip(self): - source_dirpath = Path( - iris.tests.get_data_path(["NetCDF", "unstructured_grid"]) - ) - file_name = "data_C4.nc" - source_filepath = str(source_dirpath / file_name) - - # # Ensure that our test reference text matches the original input file - # self.assertCDL(source_filepath) - - # Load the data (one cube) and save it out to a temporary file of the same name - with PARSE_UGRID_ON_LOAD.context(): - cube = iris.load_cube(source_filepath) - - print(cube) - target_filepath = str(self.temp_dir / file_name) - iris.save(cube, target_filepath) - - from subprocess import check_output - - print("") - print("ORIGINAL:") - text = check_output("ncdump -h " + source_filepath, shell=True) - print(text.decode()) - print("") - print("RE-SAVED:") - text = check_output("ncdump -h " + target_filepath, shell=True) - print(text.decode()) - # # Ensure that the saved result is identical - # self.assertCDL(target_filepath) - - # Now try loading BACK... - with PARSE_UGRID_ON_LOAD.context(): - cube_reload = iris.load_cube(target_filepath) - print("") - print("OUTPUT-RE-LOADED:") - print(cube_reload) - - mesh_orig = cube.mesh - mesh_reload = cube_reload.mesh - for propname in dir(mesh_orig): - if not propname.startswith("_"): - prop_orig = getattr(mesh_orig, propname) - if not callable(prop_orig): - prop_reload = getattr(mesh_reload, propname) - self.assertEqual( - (propname, prop_reload), (propname, prop_orig) - ) + def test_example_result_cdls(self): + # Snapshot the result of saving the example cases. + for ex_name, filepath in self.example_names_paths.items(): + target_ncfile_path = str(self.temp_dir / f"{ex_name}.nc") + # Create a netcdf file from the test CDL. + check_call( + f"ncgen {filepath} -k4 -o {target_ncfile_path}", shell=True + ) + # Fill in blank data-variables. + _add_standard_data(target_ncfile_path) + # Load as Iris data + with PARSE_UGRID_ON_LOAD.context(): + cubes = iris.load(target_ncfile_path) + # Re-save, to check the save behaviour. + resave_ncfile_path = str(self.temp_dir / f"{ex_name}_resaved.nc") + iris.save(cubes, resave_ncfile_path) + # Check the output against a CDL snapshot. + refdir_relpath = ( + "integration/experimental/ugrid_save/TestBasicSave/" + ) + reffile_name = str(Path(filepath).name).replace(".nc", ".cdl") + reffile_path = refdir_relpath + reffile_name + self.assertCDL(resave_ncfile_path, reference_filename=reffile_path) + + def test_example_roundtrips(self): + # Check that save-and-loadback leaves Iris data unchanged, + # for data derived from each UGRID example CDL. + # Snapshot the result of saving the example cases. + for ex_name, filepath in self.example_names_paths.items(): + print(f"Roundtrip checking : {ex_name}") + if "ex4" in ex_name: + # Skip this one now : still causing some problems ... + continue + target_ncfile_path = str(self.temp_dir / f"{ex_name}.nc") + # Create a netcdf file from the test CDL. + check_call( + f"ncgen {filepath} -k4 -o {target_ncfile_path}", shell=True + ) + # Fill in blank data-variables. + _add_standard_data(target_ncfile_path) + # Load the original as Iris data + with PARSE_UGRID_ON_LOAD.context(): + orig_cubes = iris.load(target_ncfile_path) + # Save-and-load-back to compare the Iris saved result. + resave_ncfile_path = str(self.temp_dir / f"{ex_name}_resaved.nc") + iris.save(orig_cubes, resave_ncfile_path) + with PARSE_UGRID_ON_LOAD.context(): + savedloaded_cubes = iris.load(resave_ncfile_path) + # This should match the original exactly + # ..EXCEPT for our inability to compare meshes. + for orig, reloaded in zip(orig_cubes, savedloaded_cubes): + for cube in (orig, reloaded): + # Remove conventions attributes, which may differ. + cube.attributes.pop("Conventions", None) + # Remove var-names, which may differ. + cube.var_name = None + # Compare the mesh contents (as we can't compare actual meshes) + self.assertEqual(orig.location, reloaded.location) + orig_mesh = orig.mesh + reloaded_mesh = reloaded.mesh + self.assertEqual( + orig_mesh.all_coords, reloaded_mesh.all_coords + ) + self.assertEqual( + orig_mesh.all_connectivities, + reloaded_mesh.all_connectivities, + ) + # Index the cubes to replace meshes with meshcoord-derived aux coords. + # This needs [:0] on the mesh dim, so do that on all dims. + keys = tuple([slice(0, None)] * orig.ndim) + orig = orig[keys] + reloaded = reloaded[keys] + # Resulting cubes, with collapsed mesh, should be IDENTICAL. + self.assertEqual(orig, reloaded) if __name__ == "__main__": diff --git a/lib/iris/tests/integration/experimental/ugrid_conventions_examples/README.txt b/lib/iris/tests/integration/experimental/ugrid_conventions_examples/README.txt new file mode 100644 index 0000000000..2f9bcfaf5e --- /dev/null +++ b/lib/iris/tests/integration/experimental/ugrid_conventions_examples/README.txt @@ -0,0 +1,46 @@ +Examples generated from CDL example sections in UGRID conventions v1.0 + ( see webpage: https://ugrid-conventions.github.io/ugrid-conventions/ ) + +CHANGES: + * all files had a data-var added, for ease of iris-roundtripping + * EX4 : + - had a couple of missing ";"s at lineends + - existing var "Mesh2_surface" is tied to the node dimension but has ':location = "face"' + - actually, ALL the formula terms should map to 'face', not nodes. + - created data-var maps (layers, faces), as reqd + +/home/h05/itpp/git/iris/iris_main/lib/iris/tests/results/ugrid_ref + $ for f in $(ls *.cdl); do n=$(echo $f | grep -o "[^.]*" | grep "ugrid"); echo $n; ncgen $n.cdl -4 -o $n.nc; done + ugrid_ex1_1d_mesh + ugrid_ex2_2d_triangular + ugrid_ex3_2d_flexible + ugrid_ex4_3d_layered + + (ncdump -h $n.nc >$n.ncdump.txt) + $ for f in $(ls *.cdl); do n=$(echo $f | grep -o "[^.]*" | grep "ugrid"); echo $n; ncdump -h $n.nc >$n.ncdump.txt; done + + (xxdiff $n.cdl $n.ncdump.txt;) + $ for f in $(ls *.cdl); do n=$(echo $f | grep -o "[^.]*" | grep "ugrid"); echo $n; xxdiff $n.cdl $n.ncdump.txt; done + +Then for compatibility testing... + (python iris_loadsave.py $n.nc) + $ for f in $(ls *.cdl); do n=$(echo $f | grep -o "[^.]*" | grep "ugrid"); echo $n; python iris_loadsave.py $n.nc; done + + + +So each of (4) examples has : + .cdl : original text from UGRID webpage + .nc : "ncgen" output = generated netcdf file + .ncdump.txt : "ncdump -h" output = re-generated CDL from nc file +(from iris_loadsave.py) + _REDUMP_cdl.txt : same as .ncdump.txt + _RESAVED.nc : from loading .nc + re-saving it + _RESAVED_REDUMP_cdl.txt : ncdump of _RESAVED.nc + +NEWSTYLE CHECKS ? +.CDL : + ==> ncgen ==> .nc + ==> load ==> ex_iris_data + ==> save ==> ex_iris_resaved + ==> load ==> ex_iris_saveload :: COMPARE with ex_iris_data + diff --git a/lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex1_1d_mesh.cdl b/lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex1_1d_mesh.cdl new file mode 100644 index 0000000000..d022fedc61 --- /dev/null +++ b/lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex1_1d_mesh.cdl @@ -0,0 +1,55 @@ +netcdf ex1_1d_mesh { +dimensions: +nMesh1_node = 5 ; // nNodes +nMesh1_edge = 4 ; // nEdges + +Two = 2; + +variables: +// Mesh topology +integer Mesh1 ; +Mesh1:cf_role = "mesh_topology" ; +Mesh1:long_name = "Topology data of 1D network" ; +Mesh1:topology_dimension = 1 ; +Mesh1:node_coordinates = "Mesh1_node_x Mesh1_node_y" ; +Mesh1:edge_node_connectivity = "Mesh1_edge_nodes" ; +Mesh1:edge_coordinates = "Mesh1_edge_x Mesh1_edge_y" ; // optional attribute +integer Mesh1_edge_nodes(nMesh1_edge, Two) ; +Mesh1_edge_nodes:cf_role = "edge_node_connectivity" ; +Mesh1_edge_nodes:long_name = "Maps every edge/link to the two nodes that it connects." ; +Mesh1_edge_nodes:start_index = 1 ; + +// Mesh node coordinates +double Mesh1_node_x(nMesh1_node) ; +Mesh1_node_x:standard_name = "longitude" ; +Mesh1_node_x:long_name = "Longitude of 1D network nodes." ; +Mesh1_node_x:units = "degrees_east" ; +double Mesh1_node_y(nMesh1_node) ; +Mesh1_node_y:standard_name = "latitude" ; +Mesh1_node_y:long_name = "Latitude of 1D network nodes." ; +Mesh1_node_y:units = "degrees_north" ; + +// Optional mesh edge coordinate variables +double Mesh1_edge_x(nMesh1_edge) ; +Mesh1_edge_x:standard_name = "longitude" ; +Mesh1_edge_x:long_name = "Characteristic longitude of 1D network edge (e.g. midpoint of the edge)." ; +Mesh1_edge_x:units = "degrees_east" ; +Mesh1_edge_x:bounds = "Mesh1_edge_xbnds" ; +double Mesh1_edge_y(nMesh1_edge) ; +Mesh1_edge_y:standard_name = "latitude" ; +Mesh1_edge_y:long_name = "Characteristic latitude of 1D network edge (e.g. midpoint of the edge)." ; +Mesh1_edge_y:units = "degrees_north" ; +Mesh1_edge_y:bounds = "Mesh1_edge_ybnds" ; +double Mesh1_edge_xbnds(nMesh1_edge,Two) ; +Mesh1_edge_xbnds:standard_name = "longitude" ; +Mesh1_edge_xbnds:long_name = "Longitude bounds of 1D network edge (i.e. begin and end longitude)." ; +Mesh1_edge_xbnds:units = "degrees_east" ; +double Mesh1_edge_ybnds(nMesh1_edge,Two) ; +Mesh1_edge_ybnds:standard_name = "latitude" ; +Mesh1_edge_ybnds:long_name = "Latitude bounds of 1D network edge (i.e. begin and end latitude)." ; +Mesh1_edge_ybnds:units = "degrees_north" ; + +float datavar(nMesh1_edge) ; + datavar:mesh = "Mesh1" ; + datavar:location = "edge" ; +} diff --git a/lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex2_2d_triangular.cdl b/lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex2_2d_triangular.cdl new file mode 100644 index 0000000000..1e4e483826 --- /dev/null +++ b/lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex2_2d_triangular.cdl @@ -0,0 +1,84 @@ +netcdf ex2_2d_triangular { +dimensions: +nMesh2_node = 4 ; // nNodes +nMesh2_edge = 5 ; // nEdges +nMesh2_face = 2 ; // nFaces + +Two = 2 ; +Three = 3 ; + +variables: +// Mesh topology +integer Mesh2 ; +Mesh2:cf_role = "mesh_topology" ; +Mesh2:long_name = "Topology data of 2D unstructured mesh" ; +Mesh2:topology_dimension = 2 ; +Mesh2:node_coordinates = "Mesh2_node_x Mesh2_node_y" ; +Mesh2:face_node_connectivity = "Mesh2_face_nodes" ; +Mesh2:face_dimension = "nMesh2_face" ; +Mesh2:edge_node_connectivity = "Mesh2_edge_nodes" ; // attribute required if variables will be defined on edges +Mesh2:edge_dimension = "nMesh2_edge" ; +Mesh2:edge_coordinates = "Mesh2_edge_x Mesh2_edge_y" ; // optional attribute (requires edge_node_connectivity) +Mesh2:face_coordinates = "Mesh2_face_x Mesh2_face_y" ; // optional attribute +Mesh2:face_edge_connectivity = "Mesh2_face_edges" ; // optional attribute (requires edge_node_connectivity) +Mesh2:face_face_connectivity = "Mesh2_face_links" ; // optional attribute +Mesh2:edge_face_connectivity = "Mesh2_edge_face_links" ; // optional attribute (requires edge_node_connectivity) +integer Mesh2_face_nodes(nMesh2_face, Three) ; +Mesh2_face_nodes:cf_role = "face_node_connectivity" ; +Mesh2_face_nodes:long_name = "Maps every triangular face to its three corner nodes." ; +Mesh2_face_nodes:start_index = 1 ; +integer Mesh2_edge_nodes(nMesh2_edge, Two) ; +Mesh2_edge_nodes:cf_role = "edge_node_connectivity" ; +Mesh2_edge_nodes:long_name = "Maps every edge to the two nodes that it connects." ; +Mesh2_edge_nodes:start_index = 1 ; + +// Optional mesh topology variables +integer Mesh2_face_edges(nMesh2_face, Three) ; +Mesh2_face_edges:cf_role = "face_edge_connectivity" ; +Mesh2_face_edges:long_name = "Maps every triangular face to its three edges." ; +Mesh2_face_edges:start_index = 1 ; +integer Mesh2_face_links(nMesh2_face, Three) ; +Mesh2_face_links:cf_role = "face_face_connectivity" ; +Mesh2_face_links:long_name = "neighbor faces for faces" ; +Mesh2_face_links:start_index = 1 ; +Mesh2_face_links:_FillValue = -999 ; +Mesh2_face_links:comment = "missing neighbor faces are indicated using _FillValue" ; +integer Mesh2_edge_face_links(nMesh2_edge, Two) ; +Mesh2_edge_face_links:cf_role = "edge_face_connectivity" ; +Mesh2_edge_face_links:long_name = "neighbor faces for edges" ; +Mesh2_edge_face_links:start_index = 1 ; +Mesh2_edge_face_links:_FillValue = -999 ; +Mesh2_edge_face_links:comment = "missing neighbor faces are indicated using _FillValue" ; + +// Mesh node coordinates +double Mesh2_node_x(nMesh2_node) ; +Mesh2_node_x:standard_name = "longitude" ; +Mesh2_node_x:long_name = "Longitude of 2D mesh nodes." ; +Mesh2_node_x:units = "degrees_east" ; +double Mesh2_node_y(nMesh2_node) ; +Mesh2_node_y:standard_name = "latitude" ; +Mesh2_node_y:long_name = "Latitude of 2D mesh nodes." ; +Mesh2_node_y:units = "degrees_north" ; + +// Optional mesh face and edge coordinate variables +double Mesh2_face_x(nMesh2_face) ; +Mesh2_face_x:standard_name = "longitude" ; +Mesh2_face_x:long_name = "Characteristics longitude of 2D mesh triangle (e.g. circumcenter coordinate)." ; +Mesh2_face_x:units = "degrees_east" ; +double Mesh2_face_y(nMesh2_face) ; +Mesh2_face_y:standard_name = "latitude" ; +Mesh2_face_y:long_name = "Characteristics latitude of 2D mesh triangle (e.g. circumcenter coordinate)." ; +Mesh2_face_y:units = "degrees_north" ; +double Mesh2_edge_x(nMesh2_edge) ; +Mesh2_edge_x:standard_name = "longitude" ; +Mesh2_edge_x:long_name = "Characteristic longitude of 2D mesh edge (e.g. midpoint of the edge)." ; +Mesh2_edge_x:units = "degrees_east" ; +double Mesh2_edge_y(nMesh2_edge) ; +Mesh2_edge_y:standard_name = "latitude" ; +Mesh2_edge_y:long_name = "Characteristic latitude of 2D mesh edge (e.g. midpoint of the edge)." ; +Mesh2_edge_y:units = "degrees_north" ; + +float datavar(nMesh2_face) ; + datavar:mesh = "Mesh2" ; + datavar:location = "face" ; +} diff --git a/lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex3_2d_flexible.cdl b/lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex3_2d_flexible.cdl new file mode 100644 index 0000000000..2fa077d152 --- /dev/null +++ b/lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex3_2d_flexible.cdl @@ -0,0 +1,99 @@ +netcdf ex3_2d_flexible { +dimensions: +nMesh2_node = 5 ; // nNodes +nMesh2_edge = 6 ; // nEdges +nMesh2_face = 2 ; // nFaces +nMaxMesh2_face_nodes = 4 ; // MaxNumNodesPerFace + +Two = 2 ; + +variables: +// Mesh topology +integer Mesh2 ; +Mesh2:cf_role = "mesh_topology" ; +Mesh2:long_name = "Topology data of 2D unstructured mesh" ; +Mesh2:topology_dimension = 2 ; +Mesh2:node_coordinates = "Mesh2_node_x Mesh2_node_y" ; +Mesh2:face_node_connectivity = "Mesh2_face_nodes" ; +Mesh2:face_dimension = "nMesh2_face" ; +Mesh2:edge_node_connectivity = "Mesh2_edge_nodes" ; // attribute required if variables will be defined on edges +Mesh2:edge_dimension = "nMesh2_edge" ; +Mesh2:edge_coordinates = "Mesh2_edge_x Mesh2_edge_y" ; // optional attribute (requires edge_node_connectivity) +Mesh2:face_coordinates = "Mesh2_face_x Mesh2_face_y" ; // optional attribute +Mesh2:face_edge_connectivity = "Mesh2_face_edges" ; // optional attribute (requires edge_node_connectivity) +Mesh2:face_face_connectivity = "Mesh2_face_links" ; // optional attribute +Mesh2:edge_face_connectivity = "Mesh2_edge_face_links" ; // optional attribute (requires edge_node_connectivity) +integer Mesh2_face_nodes(nMesh2_face, nMaxMesh2_face_nodes) ; +Mesh2_face_nodes:cf_role = "face_node_connectivity" ; +Mesh2_face_nodes:long_name = "Maps every face to its corner nodes." ; +Mesh2_face_nodes:_FillValue = 999999 ; +Mesh2_face_nodes:start_index = 1 ; +integer Mesh2_edge_nodes(nMesh2_edge, Two) ; +Mesh2_edge_nodes:cf_role = "edge_node_connectivity" ; +Mesh2_edge_nodes:long_name = "Maps every edge to the two nodes that it connects." ; +Mesh2_edge_nodes:start_index = 1 ; + +// Optional mesh topology variables +integer Mesh2_face_edges(nMesh2_face, nMaxMesh2_face_nodes) ; +Mesh2_face_edges:cf_role = "face_edge_connectivity" ; +Mesh2_face_edges:long_name = "Maps every face to its edges." ; +Mesh2_face_edges:_FillValue = 999999 ; +Mesh2_face_edges:start_index = 1 ; +integer Mesh2_face_links(nMesh2_face, nMaxMesh2_face_nodes) ; +Mesh2_face_links:cf_role = "face_face_connectivity" ; +Mesh2_face_links:long_name = "neighbor faces for faces" ; +Mesh2_face_links:start_index = 1 ; +Mesh2_face_links:_FillValue = -999 ; +Mesh2_face_links:comment = "missing edges as well as missing neighbor faces are indicated using _FillValue" ; +integer Mesh2_edge_face_links(nMesh2_edge, Two) ; +Mesh2_edge_face_links:cf_role = "edge_face_connectivity" ; +Mesh2_edge_face_links:long_name = "neighbor faces for edges" ; +Mesh2_edge_face_links:start_index = 1 ; +Mesh2_edge_face_links:_FillValue = -999 ; +Mesh2_edge_face_links:comment = "missing neighbor faces are indicated using _FillValue" ; + +// Mesh node coordinates +double Mesh2_node_x(nMesh2_node) ; +Mesh2_node_x:standard_name = "longitude" ; +Mesh2_node_x:long_name = "Longitude of 2D mesh nodes." ; +Mesh2_node_x:units = "degrees_east" ; +double Mesh2_node_y(nMesh2_node) ; +Mesh2_node_y:standard_name = "latitude" ; +Mesh2_node_y:long_name = "Latitude of 2D mesh nodes." ; +Mesh2_node_y:units = "degrees_north" ; + +// Optional mesh face and edge coordinate variables +double Mesh2_face_x(nMesh2_face) ; +Mesh2_face_x:standard_name = "longitude" ; +Mesh2_face_x:long_name = "Characteristics longitude of 2D mesh face." ; +Mesh2_face_x:units = "degrees_east" ; +Mesh2_face_x:bounds = "Mesh2_face_xbnds" ; +double Mesh2_face_y(nMesh2_face) ; +Mesh2_face_y:standard_name = "latitude" ; +Mesh2_face_y:long_name = "Characteristics latitude of 2D mesh face." ; +Mesh2_face_y:units = "degrees_north" ; +Mesh2_face_y:bounds = "Mesh2_face_ybnds" ; +double Mesh2_face_xbnds(nMesh2_face,nMaxMesh2_face_nodes) ; +Mesh2_face_xbnds:standard_name = "longitude" ; +Mesh2_face_xbnds:long_name = "Longitude bounds of 2D mesh face (i.e. corner coordinates)." ; +Mesh2_face_xbnds:units = "degrees_east" ; +Mesh2_face_xbnds:_FillValue = 9.9692099683868690E36; +double Mesh2_face_ybnds(nMesh2_face,nMaxMesh2_face_nodes) ; +Mesh2_face_ybnds:standard_name = "latitude" ; +Mesh2_face_ybnds:long_name = "Latitude bounds of 2D mesh face (i.e. corner coordinates)." ; +Mesh2_face_ybnds:units = "degrees_north" ; +Mesh2_face_ybnds:_FillValue = 9.9692099683868690E36; +double Mesh2_edge_x(nMesh2_edge) ; +Mesh2_edge_x:standard_name = "longitude" ; +Mesh2_edge_x:long_name = "Characteristic longitude of 2D mesh edge (e.g. midpoint of the edge)." ; +Mesh2_edge_x:units = "degrees_east" ; +double Mesh2_edge_y(nMesh2_edge) ; +Mesh2_edge_y:standard_name = "latitude" ; +Mesh2_edge_y:long_name = "Characteristic latitude of 2D mesh edge (e.g. midpoint of the edge)." ; +Mesh2_edge_y:units = "degrees_north" ; +// bounds variables for edges skipped + +float datavar(nMesh2_face) ; + datavar:mesh = "Mesh2" ; + datavar:location = "face" ; +} diff --git a/lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex4_3d_layered.cdl b/lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex4_3d_layered.cdl new file mode 100644 index 0000000000..5b3592c5d8 --- /dev/null +++ b/lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex4_3d_layered.cdl @@ -0,0 +1,124 @@ +netcdf ex4_3d_layered { +dimensions: +nMesh2_node = 6 ; // nNodes +nMesh2_edge = 7 ; // nEdges +nMesh2_face = 2 ; // nFaces +nMaxMesh2_face_nodes = 4 ; // MaxNumNodesPerFace +Mesh2_layers = 10 ; + +Two = 2 ; + +variables: +// Mesh topology +integer Mesh2 ; +Mesh2:cf_role = "mesh_topology" ; +Mesh2:long_name = "Topology data of 2D unstructured mesh" ; +Mesh2:topology_dimension = 2 ; +Mesh2:node_coordinates = "Mesh2_node_x Mesh2_node_y" ; +Mesh2:face_node_connectivity = "Mesh2_face_nodes" ; +Mesh2:face_dimension = "nMesh2_face" ; +Mesh2:edge_node_connectivity = "Mesh2_edge_nodes" ; // attribute required if variables will be defined on edges +Mesh2:edge_dimension = "nMesh2_edge" ; +Mesh2:edge_coordinates = "Mesh2_edge_x Mesh2_edge_y" ; // optional attribute (requires edge_node_connectivity) +Mesh2:face_coordinates = "Mesh2_face_x Mesh2_face_y" ; // optional attribute +Mesh2:face_edge_connectivity = "Mesh2_face_edges" ; // optional attribute (requires edge_node_connectivity) +Mesh2:face_face_connectivity = "Mesh2_face_links" ; // optional attribute +Mesh2:edge_face_connectivity = "Mesh2_edge_face_links" ; // optional attribute (requires edge_node_connectivity) +integer Mesh2_face_nodes(nMesh2_face, nMaxMesh2_face_nodes) ; +Mesh2_face_nodes:cf_role = "face_node_connectivity" ; +Mesh2_face_nodes:long_name = "Maps every face to its corner nodes." ; +Mesh2_face_nodes:_FillValue = 999999 ; +Mesh2_face_nodes:start_index = 1 ; +integer Mesh2_edge_nodes(nMesh2_edge, Two) ; +Mesh2_edge_nodes:cf_role = "edge_node_connectivity" ; +Mesh2_edge_nodes:long_name = "Maps every edge to the two nodes that it connects." ; +Mesh2_edge_nodes:start_index = 1 ; + +// Optional mesh topology variables +integer Mesh2_face_edges(nMesh2_face, nMaxMesh2_face_nodes) ; +Mesh2_face_edges:cf_role = "face_edge_connectivity" ; +Mesh2_face_edges:long_name = "Maps every face to its edges." ; +Mesh2_face_edges:_FillValue = 999999 ; +Mesh2_face_edges:start_index = 1 ; +integer Mesh2_face_links(nMesh2_face, nMaxMesh2_face_nodes) ; +Mesh2_face_links:cf_role = "face_face_connectivity" ; +Mesh2_face_links:long_name = "neighbor faces for faces" ; +Mesh2_face_links:start_index = 1 ; +Mesh2_face_links:_FillValue = -999 ; +Mesh2_face_links:comment = "missing edges as well as missing neighbor faces are indicated using _FillValue" ; +integer Mesh2_edge_face_links(nMesh2_edge, Two) ; +Mesh2_edge_face_links:cf_role = "edge_face_connectivity" ; +Mesh2_edge_face_links:long_name = "neighbor faces for edges" ; +Mesh2_edge_face_links:start_index = 1 ; +Mesh2_edge_face_links:_FillValue = -999 ; +Mesh2_edge_face_links:comment = "missing neighbor faces are indicated using _FillValue" ; + +// Mesh node coordinates +double Mesh2_node_x(nMesh2_node) ; +Mesh2_node_x:standard_name = "longitude" ; +Mesh2_node_x:long_name = "Longitude of 2D mesh nodes." ; +Mesh2_node_x:units = "degrees_east" ; +double Mesh2_node_y(nMesh2_node) ; +Mesh2_node_y:standard_name = "latitude" ; +Mesh2_node_y:long_name = "Latitude of 2D mesh nodes." ; +Mesh2_node_y:units = "degrees_north" ; + +// Optional mesh face and edge coordinate variables +double Mesh2_face_x(nMesh2_face) ; +Mesh2_face_x:standard_name = "longitude" ; +Mesh2_face_x:long_name = "Characteristics longitude of 2D mesh face." ; +Mesh2_face_x:units = "degrees_east" ; +Mesh2_face_x:bounds = "Mesh2_face_xbnds" ; +double Mesh2_face_y(nMesh2_face) ; +Mesh2_face_y:standard_name = "latitude" ; +Mesh2_face_y:long_name = "Characteristics latitude of 2D mesh face." ; +Mesh2_face_y:units = "degrees_north" ; +Mesh2_face_y:bounds = "Mesh2_face_ybnds" ; +double Mesh2_face_xbnds(nMesh2_face,nMaxMesh2_face_nodes) ; +Mesh2_face_xbnds:standard_name = "longitude" ; +Mesh2_face_xbnds:long_name = "Longitude bounds of 2D mesh face (i.e. corner coordinates)." ; +Mesh2_face_xbnds:units = "degrees_east" ; +Mesh2_face_xbnds:_FillValue = 9.9692099683868690E36; +double Mesh2_face_ybnds(nMesh2_face,nMaxMesh2_face_nodes) ; +Mesh2_face_ybnds:standard_name = "latitude" ; +Mesh2_face_ybnds:long_name = "Latitude bounds of 2D mesh face (i.e. corner coordinates)." ; +Mesh2_face_ybnds:units = "degrees_north" ; +Mesh2_face_ybnds:_FillValue = 9.9692099683868690E36; +double Mesh2_edge_x(nMesh2_edge) ; +Mesh2_edge_x:standard_name = "longitude" ; +Mesh2_edge_x:long_name = "Characteristic longitude of 2D mesh edge (e.g. midpoint of the edge)." ; +Mesh2_edge_x:units = "degrees_east" ; +double Mesh2_edge_y(nMesh2_edge) ; +Mesh2_edge_y:standard_name = "latitude" ; +Mesh2_edge_y:long_name = "Characteristic latitude of 2D mesh edge (e.g. midpoint of the edge)." ; +Mesh2_edge_y:units = "degrees_north" ; +// bounds variables for edges skipped + +// Vertical coordinate +double Mesh2_layers(Mesh2_layers) ; +Mesh2_layers:standard_name = "ocean_sigma_coordinate" ; +Mesh2_layers:long_name = "sigma at layer midpoints" ; +Mesh2_layers:positive = "up" ; +Mesh2_layers:formula_terms = "sigma: Mesh2_layers eta: Mesh2_surface depth: Mesh2_depth" ; +double Mesh2_depth(nMesh2_face) ; +Mesh2_depth:standard_name = "sea_floor_depth_below_geoid" ; +Mesh2_depth:units = "m" ; +Mesh2_depth:positive = "down" ; +Mesh2_depth:mesh = "Mesh2" ; +Mesh2_depth:location = "face" ; +Mesh2_depth:coordinates = "Mesh2_node_x Mesh2_node_y" ; +double Mesh2_surface(nMesh2_face) ; +Mesh2_surface:standard_name = "sea_surface_height_above_geoid" ; +Mesh2_surface:units = "m" ; +Mesh2_surface:mesh = "Mesh2" ; +Mesh2_surface:location = "face" ; +Mesh2_surface:coordinates = "Mesh2_face_x Mesh2_face_y" ; + +float datavar(Mesh2_layers, nMesh2_face) ; + datavar:mesh = "Mesh2" ; + datavar:location = "face" ; + +data: +Mesh2_layers = 0., 1., 2., 3., 4., 5., 6., 7., 8., 9. ; + +} diff --git a/lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex1_1d_mesh.cdl b/lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex1_1d_mesh.cdl new file mode 100644 index 0000000000..517991a17a --- /dev/null +++ b/lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex1_1d_mesh.cdl @@ -0,0 +1,39 @@ +dimensions: + Mesh1_edge_N_nodes = 2 ; + nMesh1_edge = 4 ; + nMesh1_node = 5 ; +variables: + int Mesh1 ; + Mesh1:cf_role = "mesh_topology" ; + Mesh1:topology_dimension = 1 ; + Mesh1:long_name = "Topology data of 1D network" ; + Mesh1:node_coordinates = "Mesh1_node_x Mesh1_node_y" ; + Mesh1:edge_coordinates = "Mesh1_edge_x Mesh1_edge_y" ; + Mesh1:edge_node_connectivity = "Mesh1_edge_nodes" ; + double Mesh1_node_x(nMesh1_node) ; + Mesh1_node_x:units = "degrees_east" ; + Mesh1_node_x:standard_name = "longitude" ; + Mesh1_node_x:long_name = "Longitude of 1D network nodes." ; + double Mesh1_node_y(nMesh1_node) ; + Mesh1_node_y:units = "degrees_north" ; + Mesh1_node_y:standard_name = "latitude" ; + Mesh1_node_y:long_name = "Latitude of 1D network nodes." ; + double Mesh1_edge_x(nMesh1_edge) ; + Mesh1_edge_x:units = "degrees_east" ; + Mesh1_edge_x:standard_name = "longitude" ; + Mesh1_edge_x:long_name = "Characteristic longitude of 1D network edge (e.g. midpoint of the edge)." ; + double Mesh1_edge_y(nMesh1_edge) ; + Mesh1_edge_y:units = "degrees_north" ; + Mesh1_edge_y:standard_name = "latitude" ; + Mesh1_edge_y:long_name = "Characteristic latitude of 1D network edge (e.g. midpoint of the edge)." ; + int Mesh1_edge_nodes(nMesh1_edge, Mesh1_edge_N_nodes) ; + Mesh1_edge_nodes:long_name = "Maps every edge/link to the two nodes that it connects." ; + Mesh1_edge_nodes:cf_role = "edge_node_connectivity" ; + Mesh1_edge_nodes:start_index = 1 ; + float datavar(nMesh1_edge) ; + datavar:mesh = "Mesh1" ; + datavar:location = "edge" ; + +// global attributes: + :Conventions = "CF-1.7" ; +} diff --git a/lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex2_2d_triangular.cdl b/lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex2_2d_triangular.cdl new file mode 100644 index 0000000000..5d01a263d6 --- /dev/null +++ b/lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex2_2d_triangular.cdl @@ -0,0 +1,75 @@ +dimensions: + Mesh2_edge_N_faces = 2 ; + Mesh2_edge_N_nodes = 2 ; + Mesh2_face_N_edges = 3 ; + Mesh2_face_N_faces = 3 ; + Mesh2_face_N_nodes = 3 ; + nMesh2_edge = 5 ; + nMesh2_face = 2 ; + nMesh2_node = 4 ; +variables: + int Mesh2 ; + Mesh2:cf_role = "mesh_topology" ; + Mesh2:topology_dimension = 2 ; + Mesh2:long_name = "Topology data of 2D unstructured mesh" ; + Mesh2:node_coordinates = "Mesh2_node_x Mesh2_node_y" ; + Mesh2:edge_coordinates = "Mesh2_edge_x Mesh2_edge_y" ; + Mesh2:face_coordinates = "Mesh2_face_x Mesh2_face_y" ; + Mesh2:face_node_connectivity = "Mesh2_face_nodes" ; + Mesh2:edge_node_connectivity = "Mesh2_edge_nodes" ; + Mesh2:face_edge_connectivity = "Mesh2_face_edges" ; + Mesh2:face_face_connectivity = "Mesh2_face_links" ; + Mesh2:edge_face_connectivity = "Mesh2_edge_face_links" ; + double Mesh2_node_x(nMesh2_node) ; + Mesh2_node_x:units = "degrees_east" ; + Mesh2_node_x:standard_name = "longitude" ; + Mesh2_node_x:long_name = "Longitude of 2D mesh nodes." ; + double Mesh2_node_y(nMesh2_node) ; + Mesh2_node_y:units = "degrees_north" ; + Mesh2_node_y:standard_name = "latitude" ; + Mesh2_node_y:long_name = "Latitude of 2D mesh nodes." ; + double Mesh2_edge_x(nMesh2_edge) ; + Mesh2_edge_x:units = "degrees_east" ; + Mesh2_edge_x:standard_name = "longitude" ; + Mesh2_edge_x:long_name = "Characteristic longitude of 2D mesh edge (e.g. midpoint of the edge)." ; + double Mesh2_edge_y(nMesh2_edge) ; + Mesh2_edge_y:units = "degrees_north" ; + Mesh2_edge_y:standard_name = "latitude" ; + Mesh2_edge_y:long_name = "Characteristic latitude of 2D mesh edge (e.g. midpoint of the edge)." ; + double Mesh2_face_x(nMesh2_face) ; + Mesh2_face_x:units = "degrees_east" ; + Mesh2_face_x:standard_name = "longitude" ; + Mesh2_face_x:long_name = "Characteristics longitude of 2D mesh triangle (e.g. circumcenter coordinate)." ; + double Mesh2_face_y(nMesh2_face) ; + Mesh2_face_y:units = "degrees_north" ; + Mesh2_face_y:standard_name = "latitude" ; + Mesh2_face_y:long_name = "Characteristics latitude of 2D mesh triangle (e.g. circumcenter coordinate)." ; + int Mesh2_face_nodes(nMesh2_face, Mesh2_face_N_nodes) ; + Mesh2_face_nodes:long_name = "Maps every triangular face to its three corner nodes." ; + Mesh2_face_nodes:cf_role = "face_node_connectivity" ; + Mesh2_face_nodes:start_index = 1 ; + int Mesh2_edge_nodes(nMesh2_edge, Mesh2_edge_N_nodes) ; + Mesh2_edge_nodes:long_name = "Maps every edge to the two nodes that it connects." ; + Mesh2_edge_nodes:cf_role = "edge_node_connectivity" ; + Mesh2_edge_nodes:start_index = 1 ; + int Mesh2_face_edges(nMesh2_face, Mesh2_face_N_edges) ; + Mesh2_face_edges:long_name = "Maps every triangular face to its three edges." ; + Mesh2_face_edges:cf_role = "face_edge_connectivity" ; + Mesh2_face_edges:start_index = 1 ; + int Mesh2_face_links(nMesh2_face, Mesh2_face_N_faces) ; + Mesh2_face_links:long_name = "neighbor faces for faces" ; + Mesh2_face_links:comment = "missing neighbor faces are indicated using _FillValue" ; + Mesh2_face_links:cf_role = "face_face_connectivity" ; + Mesh2_face_links:start_index = 1 ; + int Mesh2_edge_face_links(nMesh2_edge, Mesh2_edge_N_faces) ; + Mesh2_edge_face_links:long_name = "neighbor faces for edges" ; + Mesh2_edge_face_links:comment = "missing neighbor faces are indicated using _FillValue" ; + Mesh2_edge_face_links:cf_role = "edge_face_connectivity" ; + Mesh2_edge_face_links:start_index = 1 ; + float datavar(nMesh2_face) ; + datavar:mesh = "Mesh2" ; + datavar:location = "face" ; + +// global attributes: + :Conventions = "CF-1.7" ; +} diff --git a/lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex3_2d_flexible.cdl b/lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex3_2d_flexible.cdl new file mode 100644 index 0000000000..355799e0d8 --- /dev/null +++ b/lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex3_2d_flexible.cdl @@ -0,0 +1,75 @@ +dimensions: + Mesh2_edge_N_faces = 2 ; + Mesh2_edge_N_nodes = 2 ; + Mesh2_face_N_edges = 4 ; + Mesh2_face_N_faces = 4 ; + Mesh2_face_N_nodes = 4 ; + nMesh2_edge = 6 ; + nMesh2_face = 2 ; + nMesh2_node = 5 ; +variables: + int Mesh2 ; + Mesh2:cf_role = "mesh_topology" ; + Mesh2:topology_dimension = 2 ; + Mesh2:long_name = "Topology data of 2D unstructured mesh" ; + Mesh2:node_coordinates = "Mesh2_node_x Mesh2_node_y" ; + Mesh2:edge_coordinates = "Mesh2_edge_x Mesh2_edge_y" ; + Mesh2:face_coordinates = "Mesh2_face_x Mesh2_face_y" ; + Mesh2:face_node_connectivity = "Mesh2_face_nodes" ; + Mesh2:edge_node_connectivity = "Mesh2_edge_nodes" ; + Mesh2:face_edge_connectivity = "Mesh2_face_edges" ; + Mesh2:face_face_connectivity = "Mesh2_face_links" ; + Mesh2:edge_face_connectivity = "Mesh2_edge_face_links" ; + double Mesh2_node_x(nMesh2_node) ; + Mesh2_node_x:units = "degrees_east" ; + Mesh2_node_x:standard_name = "longitude" ; + Mesh2_node_x:long_name = "Longitude of 2D mesh nodes." ; + double Mesh2_node_y(nMesh2_node) ; + Mesh2_node_y:units = "degrees_north" ; + Mesh2_node_y:standard_name = "latitude" ; + Mesh2_node_y:long_name = "Latitude of 2D mesh nodes." ; + double Mesh2_edge_x(nMesh2_edge) ; + Mesh2_edge_x:units = "degrees_east" ; + Mesh2_edge_x:standard_name = "longitude" ; + Mesh2_edge_x:long_name = "Characteristic longitude of 2D mesh edge (e.g. midpoint of the edge)." ; + double Mesh2_edge_y(nMesh2_edge) ; + Mesh2_edge_y:units = "degrees_north" ; + Mesh2_edge_y:standard_name = "latitude" ; + Mesh2_edge_y:long_name = "Characteristic latitude of 2D mesh edge (e.g. midpoint of the edge)." ; + double Mesh2_face_x(nMesh2_face) ; + Mesh2_face_x:units = "degrees_east" ; + Mesh2_face_x:standard_name = "longitude" ; + Mesh2_face_x:long_name = "Characteristics longitude of 2D mesh face." ; + double Mesh2_face_y(nMesh2_face) ; + Mesh2_face_y:units = "degrees_north" ; + Mesh2_face_y:standard_name = "latitude" ; + Mesh2_face_y:long_name = "Characteristics latitude of 2D mesh face." ; + int Mesh2_face_nodes(nMesh2_face, Mesh2_face_N_nodes) ; + Mesh2_face_nodes:long_name = "Maps every face to its corner nodes." ; + Mesh2_face_nodes:cf_role = "face_node_connectivity" ; + Mesh2_face_nodes:start_index = 1 ; + int Mesh2_edge_nodes(nMesh2_edge, Mesh2_edge_N_nodes) ; + Mesh2_edge_nodes:long_name = "Maps every edge to the two nodes that it connects." ; + Mesh2_edge_nodes:cf_role = "edge_node_connectivity" ; + Mesh2_edge_nodes:start_index = 1 ; + int Mesh2_face_edges(nMesh2_face, Mesh2_face_N_edges) ; + Mesh2_face_edges:long_name = "Maps every face to its edges." ; + Mesh2_face_edges:cf_role = "face_edge_connectivity" ; + Mesh2_face_edges:start_index = 1 ; + int Mesh2_face_links(nMesh2_face, Mesh2_face_N_faces) ; + Mesh2_face_links:long_name = "neighbor faces for faces" ; + Mesh2_face_links:comment = "missing edges as well as missing neighbor faces are indicated using _FillValue" ; + Mesh2_face_links:cf_role = "face_face_connectivity" ; + Mesh2_face_links:start_index = 1 ; + int Mesh2_edge_face_links(nMesh2_edge, Mesh2_edge_N_faces) ; + Mesh2_edge_face_links:long_name = "neighbor faces for edges" ; + Mesh2_edge_face_links:comment = "missing neighbor faces are indicated using _FillValue" ; + Mesh2_edge_face_links:cf_role = "edge_face_connectivity" ; + Mesh2_edge_face_links:start_index = 1 ; + float datavar(nMesh2_face) ; + datavar:mesh = "Mesh2" ; + datavar:location = "face" ; + +// global attributes: + :Conventions = "CF-1.7" ; +} diff --git a/lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex4_3d_layered.cdl b/lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex4_3d_layered.cdl new file mode 100644 index 0000000000..05517f340f --- /dev/null +++ b/lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex4_3d_layered.cdl @@ -0,0 +1,106 @@ +dimensions: + Mesh2_edge_N_faces = 2 ; + Mesh2_edge_N_nodes = 2 ; + Mesh2_face_N_edges = 4 ; + Mesh2_face_N_faces = 4 ; + Mesh2_face_N_nodes = 4 ; + Mesh2_layers = 10 ; + nMesh2_edge = 7 ; + nMesh2_face = 2 ; + nMesh2_node = 6 ; +variables: + int Mesh2 ; + Mesh2:cf_role = "mesh_topology" ; + Mesh2:topology_dimension = 2 ; + Mesh2:long_name = "Topology data of 2D unstructured mesh" ; + Mesh2:node_coordinates = "Mesh2_node_x Mesh2_node_y" ; + Mesh2:edge_coordinates = "Mesh2_edge_x Mesh2_edge_y" ; + Mesh2:face_coordinates = "Mesh2_face_x Mesh2_face_y" ; + Mesh2:face_node_connectivity = "Mesh2_face_nodes" ; + Mesh2:edge_node_connectivity = "Mesh2_edge_nodes" ; + Mesh2:face_edge_connectivity = "Mesh2_face_edges" ; + Mesh2:face_face_connectivity = "Mesh2_face_links" ; + Mesh2:edge_face_connectivity = "Mesh2_edge_face_links" ; + double Mesh2_node_x(nMesh2_node) ; + Mesh2_node_x:units = "degrees_east" ; + Mesh2_node_x:standard_name = "longitude" ; + Mesh2_node_x:long_name = "Longitude of 2D mesh nodes." ; + double Mesh2_node_y(nMesh2_node) ; + Mesh2_node_y:units = "degrees_north" ; + Mesh2_node_y:standard_name = "latitude" ; + Mesh2_node_y:long_name = "Latitude of 2D mesh nodes." ; + double Mesh2_edge_x(nMesh2_edge) ; + Mesh2_edge_x:units = "degrees_east" ; + Mesh2_edge_x:standard_name = "longitude" ; + Mesh2_edge_x:long_name = "Characteristic longitude of 2D mesh edge (e.g. midpoint of the edge)." ; + double Mesh2_edge_y(nMesh2_edge) ; + Mesh2_edge_y:units = "degrees_north" ; + Mesh2_edge_y:standard_name = "latitude" ; + Mesh2_edge_y:long_name = "Characteristic latitude of 2D mesh edge (e.g. midpoint of the edge)." ; + double Mesh2_face_x(nMesh2_face) ; + Mesh2_face_x:units = "degrees_east" ; + Mesh2_face_x:standard_name = "longitude" ; + Mesh2_face_x:long_name = "Characteristics longitude of 2D mesh face." ; + double Mesh2_face_y(nMesh2_face) ; + Mesh2_face_y:units = "degrees_north" ; + Mesh2_face_y:standard_name = "latitude" ; + Mesh2_face_y:long_name = "Characteristics latitude of 2D mesh face." ; + int Mesh2_face_nodes(nMesh2_face, Mesh2_face_N_nodes) ; + Mesh2_face_nodes:long_name = "Maps every face to its corner nodes." ; + Mesh2_face_nodes:cf_role = "face_node_connectivity" ; + Mesh2_face_nodes:start_index = 1 ; + int Mesh2_edge_nodes(nMesh2_edge, Mesh2_edge_N_nodes) ; + Mesh2_edge_nodes:long_name = "Maps every edge to the two nodes that it connects." ; + Mesh2_edge_nodes:cf_role = "edge_node_connectivity" ; + Mesh2_edge_nodes:start_index = 1 ; + int Mesh2_face_edges(nMesh2_face, Mesh2_face_N_edges) ; + Mesh2_face_edges:long_name = "Maps every face to its edges." ; + Mesh2_face_edges:cf_role = "face_edge_connectivity" ; + Mesh2_face_edges:start_index = 1 ; + int Mesh2_face_links(nMesh2_face, Mesh2_face_N_faces) ; + Mesh2_face_links:long_name = "neighbor faces for faces" ; + Mesh2_face_links:comment = "missing edges as well as missing neighbor faces are indicated using _FillValue" ; + Mesh2_face_links:cf_role = "face_face_connectivity" ; + Mesh2_face_links:start_index = 1 ; + int Mesh2_edge_face_links(nMesh2_edge, Mesh2_edge_N_faces) ; + Mesh2_edge_face_links:long_name = "neighbor faces for edges" ; + Mesh2_edge_face_links:comment = "missing neighbor faces are indicated using _FillValue" ; + Mesh2_edge_face_links:cf_role = "edge_face_connectivity" ; + Mesh2_edge_face_links:start_index = 1 ; + float datavar(Mesh2_layers, nMesh2_face) ; + datavar:mesh = "Mesh2" ; + datavar:location = "face" ; + datavar:coordinates = "Mesh2_depth Mesh2_surface" ; + double Mesh2_layers(Mesh2_layers) ; + Mesh2_layers:axis = "Z" ; + Mesh2_layers:units = "1" ; + Mesh2_layers:standard_name = "ocean_sigma_coordinate" ; + Mesh2_layers:long_name = "sigma at layer midpoints" ; + Mesh2_layers:positive = "up" ; + Mesh2_layers:formula_terms = "sigma: Mesh2_layers eta: Mesh2_surface depth: Mesh2_depth" ; + double Mesh2_depth(nMesh2_face) ; + Mesh2_depth:units = "m" ; + Mesh2_depth:standard_name = "sea_floor_depth_below_geoid" ; + Mesh2_depth:location = "face" ; + Mesh2_depth:mesh = "Mesh2" ; + Mesh2_depth:positive = "down" ; + double Mesh2_surface(nMesh2_face) ; + Mesh2_surface:units = "m" ; + Mesh2_surface:standard_name = "sea_surface_height_above_geoid" ; + Mesh2_surface:location = "face" ; + Mesh2_surface:mesh = "Mesh2" ; + double Mesh2_depth_0(nMesh2_face) ; + Mesh2_depth_0:standard_name = "sea_floor_depth_below_geoid" ; + Mesh2_depth_0:units = "m" ; + Mesh2_depth_0:positive = "down" ; + Mesh2_depth_0:mesh = "Mesh2" ; + Mesh2_depth_0:location = "face" ; + double Mesh2_surface_0(nMesh2_face) ; + Mesh2_surface_0:standard_name = "sea_surface_height_above_geoid" ; + Mesh2_surface_0:units = "m" ; + Mesh2_surface_0:mesh = "Mesh2" ; + Mesh2_surface_0:location = "face" ; + +// global attributes: + :Conventions = "CF-1.7" ; +} diff --git a/lib/iris/tests/stock/netcdf.py b/lib/iris/tests/stock/netcdf.py index 78dc08eafd..30fd0facce 100644 --- a/lib/iris/tests/stock/netcdf.py +++ b/lib/iris/tests/stock/netcdf.py @@ -62,6 +62,8 @@ def _add_standard_data(nc_path, unlimited_dim_size=0): ] # Data addition dependent on this assumption: assert len(unlimited_dim_names) < 2 + if len(unlimited_dim_names) == 0: + unlimited_dim_names = ["*unused*"] # Fill variables data with placeholder numbers. for var in ds.variables.values(): @@ -72,11 +74,13 @@ def _add_standard_data(nc_path, unlimited_dim_size=0): unlimited_dim_size if dim == unlimited_dim_names[0] else size for dim, size in zip(dims, shape) ] - data = np.zeros(shape, dtype=var.dtype) + data = np.ones(shape, dtype=var.dtype) # Do not use zero if len(var.dimensions) == 1 and var.dimensions[0] == var.name: # Fill the var with ascending values (not all zeroes), # so it can be a dim-coord. - data = np.arange(data.size, dtype=data.dtype).reshape(data.shape) + data = np.arange(1, data.size + 1, dtype=data.dtype).reshape( + data.shape + ) var[:] = data ds.close() From f5c471ca4531cc573969512219494a45f2d7d124 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 2 Sep 2021 20:59:24 +0100 Subject: [PATCH 07/26] Remove bad comment. --- lib/iris/tests/integration/experimental/test_ugrid_save.py | 1 - 1 file changed, 1 deletion(-) diff --git a/lib/iris/tests/integration/experimental/test_ugrid_save.py b/lib/iris/tests/integration/experimental/test_ugrid_save.py index 969b5d7557..236ba006b0 100644 --- a/lib/iris/tests/integration/experimental/test_ugrid_save.py +++ b/lib/iris/tests/integration/experimental/test_ugrid_save.py @@ -71,7 +71,6 @@ def test_example_result_cdls(self): def test_example_roundtrips(self): # Check that save-and-loadback leaves Iris data unchanged, # for data derived from each UGRID example CDL. - # Snapshot the result of saving the example cases. for ex_name, filepath in self.example_names_paths.items(): print(f"Roundtrip checking : {ex_name}") if "ex4" in ex_name: From 433efb06eeb84b87a2304b7a9bbdbf8398ca602e Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 2 Sep 2021 21:32:13 +0100 Subject: [PATCH 08/26] Start on specific tests. --- .../fileformats/netcdf/test_Saver__ugrid.py | 213 ++++++++++++++++++ 1 file changed, 213 insertions(+) create mode 100644 lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py diff --git a/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py b/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py new file mode 100644 index 0000000000..f08e1218a2 --- /dev/null +++ b/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py @@ -0,0 +1,213 @@ +# Copyright Iris contributors +# +# This file is part of Iris and is released under the LGPL license. +# See COPYING and COPYING.LESSER in the root of the repository for full +# licensing details. +"""Unit tests for the `iris.fileformats.netcdf.Saver` class.""" + +# Import iris.tests first so that some things can be initialised before +# importing anything else. +import iris.tests as tests # isort:skip + +from pathlib import Path +import shutil +import tempfile + +import numpy as np + +from iris import save +from iris.coords import AuxCoord +from iris.cube import Cube +from iris.experimental.ugrid import Connectivity, Mesh + +XY_LOCS = ("x", "y") +XY_NAMES = ("longitude", "latitude") + + +def build_mesh( + n_nodes=2, + n_faces=0, + n_edges=0, + nodecoord_xyargs=None, + edgecoord_xyargs=None, + facecoord_xyargs=None, + conn_role_kwargs={}, # role: kwargs + mesh_kwargs=None, +): + """ + Make a test mesh. + + Mesh has faces edges, face-coords and edge-coords, numbers of which can be controlled. + + """ + + def applyargs(coord, kwargs): + if kwargs: + for key, val in kwargs.items(): + # kwargs is a dict + setattr(coord, key, val) + + def apply_xyargs(coords, xyargs): + if xyargs: + for coord, kwargs in zip(coords, xyargs): + # coords and xyargs both iterables : implicitly=(x,y) + applyargs(coord, kwargs) + + # NB when creating coords, supply axis to make Mesh.to_AuxCoords work + node_coords = [ + AuxCoord(np.arange(n_nodes), standard_name=name) + for loc, name in zip(XY_LOCS, XY_NAMES) + ] + apply_xyargs(node_coords, nodecoord_xyargs) + + connectivities = {} + edge_coords = [] + face_coords = [] + topology_dimension = 0 + if n_edges: + topology_dimension = 1 + connectivities["edge_node_connectivity"] = Connectivity( + np.zeros((n_edges, 2), np.int32), cf_role="edge_node_connectivity" + ) + edge_coords = [ + AuxCoord(np.arange(n_edges), standard_name=name) + for loc, name in zip(XY_LOCS, XY_NAMES) + ] + apply_xyargs(edge_coords, edgecoord_xyargs) + + if n_faces: + topology_dimension = 2 + connectivities["face_node_connectivity"] = Connectivity( + np.zeros((n_faces, 4), np.int32), cf_role="face_node_connectivity" + ) + face_coords = [ + AuxCoord(np.arange(n_faces), standard_name=name) + for loc, name in zip(XY_LOCS, XY_NAMES) + ] + apply_xyargs(face_coords, facecoord_xyargs) + + mesh_dims = {"node": n_nodes, "edge": n_edges, "face": n_faces} + + for role, kwargs in conn_role_kwargs.items(): + if role in connectivities: + conn = connectivities[role] + else: + loc_from, loc_to, _ = role.split("_") + dims = [mesh_dims[loc] for loc in (loc_from, loc_to)] + conn = Connectivity(np.zeros(dims, dtype=np.int32), cf_role=role) + connectivities[role] = conn + applyargs(conn, kwargs) + + mesh = Mesh( + topology_dimension=topology_dimension, + node_coords_and_axes=zip(node_coords, XY_LOCS), + edge_coords_and_axes=zip(edge_coords, XY_LOCS), + face_coords_and_axes=zip(face_coords, XY_LOCS), + connectivities=connectivities.values(), + ) + applyargs(mesh, mesh_kwargs) + mesh._mesh_dims = mesh_dims + + return mesh + + +def make_mesh(basic=True, **kwargs): + if basic: + # Use some helpful non-minimal settings as our 'basic' mesh. + use_kwargs = dict( + n_nodes=5, + n_faces=2, + nodecoord_xyargs=tuple( + dict(var_name=f"node_{loc}") for loc in XY_LOCS + ), + facecoord_xyargs=tuple( + dict(var_name=f"face_{loc}") for loc in XY_LOCS + ), + mesh_kwargs=dict( + var_name="Mesh2d", + node_dimension="Mesh2d_nodes", + face_dimension="Mesh2d_faces", + ), + ) + use_kwargs.update(kwargs) + else: + use_kwargs = kwargs + + mesh = build_mesh(**use_kwargs) + return mesh + + +# Make simple "standard" test meshes for multiple uses +_DEFAULT_MESH = make_mesh() +# NB the 'minimal' mesh might have just nodes, but Mesh doesn't (yet) support it +_MINIMAL_MESH = make_mesh(basic=False, n_nodes=3, n_edges=2, n_faces=0) + + +def make_cube(mesh=_DEFAULT_MESH, location="face", **kwargs): + dim = mesh._mesh_dims[location] + cube = Cube(np.zeros(dim, np.float32)) + for meshco in mesh.to_MeshCoords(location): + cube.add_aux_coord(meshco, (0,)) + for key, val in kwargs.items(): + setattr(cube, key, val) + return cube + + +from subprocess import check_output + + +class TestSaveUgrid__cube(tests.IrisTest): + @classmethod + def setUpClass(cls): + cls.temp_dir = Path(tempfile.mkdtemp()) + cls.tempfile_path = str(cls.temp_dir / "tmp.nc") + + @classmethod + def tearDownClass(cls): + shutil.rmtree(cls.temp_dir) + + def check_save(self, data): + save(data, self.tempfile_path) + text = check_output(f"ncdump -h {self.tempfile_path}", shell=True) + text = text.decode() + print(text) + + def test_minimal_mesh_cdl(self): + data = make_cube(mesh=_MINIMAL_MESH, location="edge") + self.check_save(data) + # self.assertCDL(self.tempfile_path) + + def test_basic_mesh_cdl(self): + data = make_cube(mesh=_DEFAULT_MESH) + self.check_save(data) + # self.assertCDL(self.tempfile_path) + + def test_multi_cubes_common_mesh(self): + cube1 = make_cube(var_name="a") + cube2 = make_cube(var_name="b") + self.check_save([cube1, cube2]) + # self.assertCDL(self.tempfile_path) + + def test_multi_cubes_identical_mesh(self): + mesh1 = make_mesh() + mesh2 = make_mesh() + cube1 = make_cube(var_name="a", mesh=mesh1) + cube2 = make_cube(var_name="b", mesh=mesh2) + self.check_save([cube1, cube2]) + # self.assertCDL(self.tempfile_path) + + def test_multi_cubes_different_mesh(self): + cube1 = make_cube(var_name="a") + cube2 = make_cube(var_name="b", mesh=make_mesh(n_faces=4)) + self.check_save([cube1, cube2]) + # self.assertCDL(self.tempfile_path) + + def test_multi_cubes_different_locations(self): + cube1 = make_cube(var_name="a", location="face") + cube2 = make_cube(var_name="b", location="node") + self.check_save([cube1, cube2]) + # self.assertCDL(self.tempfile_path) + + +if __name__ == "__main__": + tests.main() From f126860e367f8446372d5ad2fd8a29e1e0f30ba8 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Fri, 3 Sep 2021 01:54:36 +0100 Subject: [PATCH 09/26] First proper result checking. --- .../fileformats/netcdf/test_Saver__ugrid.py | 132 +++++++++++++++++- 1 file changed, 126 insertions(+), 6 deletions(-) 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 f08e1218a2..ccd15cdb65 100644 --- a/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py +++ b/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py @@ -11,6 +11,7 @@ from pathlib import Path import shutil +from subprocess import check_output import tempfile import numpy as np @@ -153,7 +154,70 @@ def make_cube(mesh=_DEFAULT_MESH, location="face", **kwargs): return cube -from subprocess import check_output +def scan_dataset(filepath): + """ + Snapshot a dataset (the key metadata). + + Returns: + a tuple (dimsdict, varsdict) + * dimsdict (dict): + A mapping of dimension-name: length. + * varsdict (dict): + A map of each variable's properties, {var_name: propsdict} + Each propsdict is {attribute-name: value} over the var's ncattrs(). + Each propsdict ALSO contains a ['_dims'] entry listing the + variable's dims. + + """ + import netCDF4 as nc + + ds = nc.Dataset(filepath) + # dims dict is {name: len} + dimsdict = {name: dim.size for name, dim in ds.dimensions.items()} + # vars dict is {name: {attr:val}} + varsdict = {} + for name, var in ds.variables.items(): + varsdict[name] = {prop: getattr(var, prop) for prop in var.ncattrs()} + varsdict[name]["_dims"] = list(var.dimensions) + ds.close() + return dimsdict, varsdict + + +def vars_w_props(varsdict, **kwargs): + """ + Subset a vars dict, {name:props}, for those where given attributes=values. + Except '="*"' means the '' attribute must only exist. + + """ + + def check_attrs_match(attrs): + result = True + for key, val in kwargs.items(): + result = key in attrs + if result and val != "*": + # val='*'' for a simple existence check + # Otherwise actual value must also match + result = attrs[key] == val + if not result: + break + return result + + varsdict = { + name: attrs + for name, attrs in varsdict.items() + if check_attrs_match(attrs) + } + return varsdict + + +def vars_w_dims(varsdict, dim_names): + """Subset a vars dict elect, selecting only those with all the given dimensions.""" + varsdict = { + name: propsdict + for name, propsdict in varsdict.items() + if all(dim in propsdict["_dims"] for dim in dim_names) + } + return varsdict class TestSaveUgrid__cube(tests.IrisTest): @@ -171,11 +235,34 @@ def check_save(self, data): text = check_output(f"ncdump -h {self.tempfile_path}", shell=True) text = text.decode() print(text) + return scan_dataset(self.tempfile_path) def test_minimal_mesh_cdl(self): data = make_cube(mesh=_MINIMAL_MESH, location="edge") - self.check_save(data) - # self.assertCDL(self.tempfile_path) + dims, vars = self.check_save(data) + + # There should be 1 mesh var. + mesh_vars = vars_w_props(vars, cf_role="mesh_topology") + self.assertEqual(1, len(mesh_vars)) + (mesh_name,) = mesh_vars.keys() + + # There should be 1 mesh-linked (data)var + data_vars = vars_w_props(vars, mesh="*") + self.assertEqual(1, len(data_vars)) + + # The mesh var should link to the mesh, at 'edges' + self.assertEqual(["unknown"], list(data_vars.keys())) + (a_props,) = data_vars.values() + self.assertEqual(mesh_name, a_props["mesh"]) + self.assertEqual("edge", a_props["location"]) + + # get name of first edge coord + edge_coord = vars[mesh_name]["edge_coordinates"].split(" ")[0] + # get edge dim = first dim of edge coord + (edge_dim,) = vars[edge_coord]["_dims"] + + # The dims of the datavar should == [edges] + self.assertEqual([edge_dim], a_props["_dims"]) def test_basic_mesh_cdl(self): data = make_cube(mesh=_DEFAULT_MESH) @@ -185,21 +272,54 @@ def test_basic_mesh_cdl(self): def test_multi_cubes_common_mesh(self): cube1 = make_cube(var_name="a") cube2 = make_cube(var_name="b") - self.check_save([cube1, cube2]) + dims, vars = self.check_save([cube1, cube2]) + # there is only 1 mesh in the file + mesh_vars = vars_w_props(vars, cf_role="mesh_topology") + self.assertEqual(len(mesh_vars), 1) + (mesh_name,) = mesh_vars.keys() + # both the main variables reference the same mesh, and 'face' location + v_a, v_b = vars["a"], vars["b"] + self.assertEqual(mesh_name, v_a["mesh"]) + self.assertEqual("face", v_a["location"]) + self.assertEqual(mesh_name, v_b["mesh"]) + self.assertEqual("face", v_b["location"]) # self.assertCDL(self.tempfile_path) - def test_multi_cubes_identical_mesh(self): + def test_multi_cubes_identical_meshes(self): mesh1 = make_mesh() mesh2 = make_mesh() cube1 = make_cube(var_name="a", mesh=mesh1) cube2 = make_cube(var_name="b", mesh=mesh2) - self.check_save([cube1, cube2]) + cube1 = make_cube(var_name="a") + cube2 = make_cube(var_name="b") + dims, vars = self.check_save([cube1, cube2]) + # there are 2 meshes in the file + mesh_vars = vars_w_props(vars, cf_role="mesh_topology") + self.assertEqual(len(mesh_vars), 2) + # there are two (data)variables with a 'mesh' property + mesh_datavars = vars_w_props(vars, mesh="*") + self.assertEqual(2, len(mesh_datavars)) + self.assertEqual(["a", "b"], sorted(mesh_datavars.keys())) + # the main variables reference the same mesh, and 'face' location + v_a, v_b = vars["a"], vars["b"] + mesh_a, mesh_b = v_a["mesh"], v_b["mesh"] + self.assertNotEqual(mesh_a, mesh_b) + self.assertEqual(sorted([mesh_a, mesh_b]), sorted(mesh_vars.keys())) # self.assertCDL(self.tempfile_path) def test_multi_cubes_different_mesh(self): cube1 = make_cube(var_name="a") cube2 = make_cube(var_name="b", mesh=make_mesh(n_faces=4)) self.check_save([cube1, cube2]) + dims, vars = self.check_save([cube1, cube2]) + # there are 2 meshes in the file + mesh_vars = vars_w_props(vars, cf_role="mesh_topology") + self.assertEqual(len(mesh_vars), 2) + # there are two (data)variables with a 'mesh' property + mesh_datavars = vars_w_props(vars, mesh="*") + self.assertEqual(2, len(mesh_datavars)) + self.assertEqual(["a", "b"], sorted(mesh_datavars.keys())) + # the main variables reference the same mesh, and 'face' location # self.assertCDL(self.tempfile_path) def test_multi_cubes_different_locations(self): From 532bb8b53c4142297291fb25ceb54b1764788601 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Fri, 3 Sep 2021 17:21:27 +0100 Subject: [PATCH 10/26] Slight tidy, drop _MINIMAL_MESH, new testing util routines. --- .../fileformats/netcdf/test_Saver__ugrid.py | 192 ++++++++++++------ 1 file changed, 125 insertions(+), 67 deletions(-) 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 ccd15cdb65..978c641269 100644 --- a/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py +++ b/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py @@ -32,7 +32,7 @@ def build_mesh( nodecoord_xyargs=None, edgecoord_xyargs=None, facecoord_xyargs=None, - conn_role_kwargs={}, # role: kwargs + conn_role_kwargs=None, # mapping {connectivity-role: connectivity-kwargs} mesh_kwargs=None, ): """ @@ -89,15 +89,18 @@ def apply_xyargs(coords, xyargs): mesh_dims = {"node": n_nodes, "edge": n_edges, "face": n_faces} - for role, kwargs in conn_role_kwargs.items(): - if role in connectivities: - conn = connectivities[role] - else: - loc_from, loc_to, _ = role.split("_") - dims = [mesh_dims[loc] for loc in (loc_from, loc_to)] - conn = Connectivity(np.zeros(dims, dtype=np.int32), cf_role=role) - connectivities[role] = conn - applyargs(conn, kwargs) + if conn_role_kwargs: + for role, kwargs in conn_role_kwargs.items(): + if role in connectivities: + conn = connectivities[role] + else: + loc_from, loc_to, _ = role.split("_") + dims = [mesh_dims[loc] for loc in (loc_from, loc_to)] + conn = Connectivity( + np.zeros(dims, dtype=np.int32), cf_role=role + ) + connectivities[role] = conn + applyargs(conn, kwargs) mesh = Mesh( topology_dimension=topology_dimension, @@ -138,10 +141,8 @@ def make_mesh(basic=True, **kwargs): return mesh -# Make simple "standard" test meshes for multiple uses +# Pre-create a simple "standard" test mesh for multiple uses _DEFAULT_MESH = make_mesh() -# NB the 'minimal' mesh might have just nodes, but Mesh doesn't (yet) support it -_MINIMAL_MESH = make_mesh(basic=False, n_nodes=3, n_edges=2, n_faces=0) def make_cube(mesh=_DEFAULT_MESH, location="face", **kwargs): @@ -156,17 +157,17 @@ def make_cube(mesh=_DEFAULT_MESH, location="face", **kwargs): def scan_dataset(filepath): """ - Snapshot a dataset (the key metadata). + Snapshot a netcdf dataset (the key metadata). Returns: - a tuple (dimsdict, varsdict) + dimsdict, varsdict * dimsdict (dict): - A mapping of dimension-name: length. + A map of dimension-name: length. * varsdict (dict): A map of each variable's properties, {var_name: propsdict} Each propsdict is {attribute-name: value} over the var's ncattrs(). Each propsdict ALSO contains a ['_dims'] entry listing the - variable's dims. + variable's dims. """ import netCDF4 as nc @@ -185,8 +186,9 @@ def scan_dataset(filepath): def vars_w_props(varsdict, **kwargs): """ - Subset a vars dict, {name:props}, for those where given attributes=values. - Except '="*"' means the '' attribute must only exist. + Subset a vars dict, {name:props}, returning only those where each + =, defined by the given keywords. + Except that '="*"' means that an attribute '' merely _exists_. """ @@ -194,10 +196,9 @@ def check_attrs_match(attrs): result = True for key, val in kwargs.items(): result = key in attrs - if result and val != "*": + if result: # val='*'' for a simple existence check - # Otherwise actual value must also match - result = attrs[key] == val + result = (val == "*") or attrs[key] == val if not result: break return result @@ -211,7 +212,7 @@ def check_attrs_match(attrs): def vars_w_dims(varsdict, dim_names): - """Subset a vars dict elect, selecting only those with all the given dimensions.""" + """Subset a vars dict, returning all those which map all the specified dims.""" varsdict = { name: propsdict for name, propsdict in varsdict.items() @@ -220,6 +221,38 @@ def vars_w_dims(varsdict, dim_names): return varsdict +def vars_meshvars(vars): + """Subset a varsdict, returning those which are mesh variables (by cf_role).""" + return vars_w_props(vars, cf_role="mesh_topology") + + +def vars_meshdim(vars, location, mesh_name=None): + """ " + Extract a dim-name for a given element location. + + Args: + * vars (varsdict): + file varsdict, as returned from 'snapshot_dataset'. + * location (string): + a mesh location : 'node' / 'edge' / 'face' + * mesh_name (string or None): + If given, identifies the mesh var. + Otherwise, find a unique mesh var (i.e. there must be exactly 1). + + Returns: + dim_name (string) + The dim-name of the mesh dim for the given location. + + """ + if mesh_name is None: + # Find "the" meshvar -- assuming there is just one. + (mesh_name,) = vars_meshvars(vars).keys() + mesh_props = vars[mesh_name] + loc_coords = mesh_props[f"{location}_coordinates"].split(" ") + (single_location_dim,) = vars[loc_coords[0]]["_dims"] + return single_location_dim + + class TestSaveUgrid__cube(tests.IrisTest): @classmethod def setUpClass(cls): @@ -237,54 +270,84 @@ def check_save(self, data): print(text) return scan_dataset(self.tempfile_path) - def test_minimal_mesh_cdl(self): - data = make_cube(mesh=_MINIMAL_MESH, location="edge") + def test_basic_mesh_cdl(self): + data = make_cube(mesh=_DEFAULT_MESH) dims, vars = self.check_save(data) - # There should be 1 mesh var. - mesh_vars = vars_w_props(vars, cf_role="mesh_topology") - self.assertEqual(1, len(mesh_vars)) - (mesh_name,) = mesh_vars.keys() + # There is exactly 1 mesh var. + ((mesh_name, mesh_props),) = vars_meshvars(vars).items() - # There should be 1 mesh-linked (data)var + # There is exactly 1 mesh-linked (data)var data_vars = vars_w_props(vars, mesh="*") - self.assertEqual(1, len(data_vars)) - - # The mesh var should link to the mesh, at 'edges' - self.assertEqual(["unknown"], list(data_vars.keys())) - (a_props,) = data_vars.values() - self.assertEqual(mesh_name, a_props["mesh"]) - self.assertEqual("edge", a_props["location"]) + ((a_name, a_props),) = data_vars.items() + + # The mesh var links to the mesh, with location 'faces' + self.assertEqual(a_name, "unknown") + self.assertEqual(a_props["mesh"], mesh_name) + self.assertEqual(a_props["location"], "face") + + # There are 2 face coords == those listed in the mesh + face_coords = mesh_props["face_coordinates"].split(" ") + self.assertEqual(len(face_coords), 2) + # get face dim (actually, from the dims of the first face coord) + face_dim = vars_meshdim(vars, "face") + # The face coords should both map that single dim. + self.assertTrue( + all(vars[co]["_dims"] == [face_dim] for co in face_coords) + ) - # get name of first edge coord - edge_coord = vars[mesh_name]["edge_coordinates"].split(" ")[0] - # get edge dim = first dim of edge coord - (edge_dim,) = vars[edge_coord]["_dims"] + # The dims of the datavar also == [] + self.assertEqual(a_props["_dims"], [face_dim]) - # The dims of the datavar should == [edges] - self.assertEqual([edge_dim], a_props["_dims"]) + # There are 2 node coordinates == those listed in the mesh. + node_coords = mesh_props["node_coordinates"].split(" ") + self.assertEqual(len(node_coords), 2) + # These should be the only ones using the 'nodes' dimension. + node_dim = vars_meshdim(vars, "node") + self.assertEqual( + sorted(node_coords), sorted(vars_w_dims(vars, [node_dim]).keys()) + ) - def test_basic_mesh_cdl(self): - data = make_cube(mesh=_DEFAULT_MESH) - self.check_save(data) - # self.assertCDL(self.tempfile_path) + # There are no edges. + self.assertEqual( + len(vars_w_props(vars, cf_role="edge_node_connectivity")), 0 + ) def test_multi_cubes_common_mesh(self): cube1 = make_cube(var_name="a") cube2 = make_cube(var_name="b") dims, vars = self.check_save([cube1, cube2]) - # there is only 1 mesh in the file - mesh_vars = vars_w_props(vars, cf_role="mesh_topology") - self.assertEqual(len(mesh_vars), 1) - (mesh_name,) = mesh_vars.keys() + # there is exactly 1 mesh in the file + ((mesh_name, _),) = vars_meshvars(vars).items() # both the main variables reference the same mesh, and 'face' location v_a, v_b = vars["a"], vars["b"] - self.assertEqual(mesh_name, v_a["mesh"]) - self.assertEqual("face", v_a["location"]) - self.assertEqual(mesh_name, v_b["mesh"]) - self.assertEqual("face", v_b["location"]) + self.assertEqual(v_a["mesh"], mesh_name) + self.assertEqual(v_a["location"], "face") + self.assertEqual(v_b["mesh"], mesh_name) + self.assertEqual(v_b["location"], "face") # self.assertCDL(self.tempfile_path) + def test_multi_cubes_different_locations(self): + cube1 = make_cube(var_name="a", location="face") + cube2 = make_cube(var_name="b", location="node") + dims, vars = self.check_save([cube1, cube2]) + + # there is exactly 1 mesh in the file + ((mesh_name, mesh_props),) = vars_meshvars(vars).items() + + # the main variables reference the same mesh at different locations + v_a, v_b = vars["a"], vars["b"] + self.assertEqual(v_a["mesh"], mesh_name) + self.assertEqual(v_a["location"], "face") + self.assertEqual(v_b["mesh"], mesh_name) + self.assertEqual(v_b["location"], "node") + + # the main variables map the face and node dimensions + face_dim = vars_meshdim(vars, "face") + node_dim = vars_meshdim(vars, "node") + self.assertEqual(v_a["_dims"], [face_dim]) + self.assertEqual(v_b["_dims"], [node_dim]) + def test_multi_cubes_identical_meshes(self): mesh1 = make_mesh() mesh2 = make_mesh() @@ -293,19 +356,20 @@ def test_multi_cubes_identical_meshes(self): cube1 = make_cube(var_name="a") cube2 = make_cube(var_name="b") dims, vars = self.check_save([cube1, cube2]) + # there are 2 meshes in the file - mesh_vars = vars_w_props(vars, cf_role="mesh_topology") + mesh_vars = vars_meshvars(vars) self.assertEqual(len(mesh_vars), 2) + # there are two (data)variables with a 'mesh' property mesh_datavars = vars_w_props(vars, mesh="*") self.assertEqual(2, len(mesh_datavars)) self.assertEqual(["a", "b"], sorted(mesh_datavars.keys())) - # the main variables reference the same mesh, and 'face' location - v_a, v_b = vars["a"], vars["b"] - mesh_a, mesh_b = v_a["mesh"], v_b["mesh"] - self.assertNotEqual(mesh_a, mesh_b) + + # the data variables reference the two separate meshes + a_props, b_props = vars["a"], vars["b"] + mesh_a, mesh_b = a_props["mesh"], b_props["mesh"] self.assertEqual(sorted([mesh_a, mesh_b]), sorted(mesh_vars.keys())) - # self.assertCDL(self.tempfile_path) def test_multi_cubes_different_mesh(self): cube1 = make_cube(var_name="a") @@ -322,12 +386,6 @@ def test_multi_cubes_different_mesh(self): # the main variables reference the same mesh, and 'face' location # self.assertCDL(self.tempfile_path) - def test_multi_cubes_different_locations(self): - cube1 = make_cube(var_name="a", location="face") - cube2 = make_cube(var_name="b", location="node") - self.check_save([cube1, cube2]) - # self.assertCDL(self.tempfile_path) - if __name__ == "__main__": tests.main() From 3d800a7a4f2c9e824a15c7190dbc6d469ffe91f9 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Mon, 6 Sep 2021 11:30:17 +0100 Subject: [PATCH 11/26] Small changes. --- lib/iris/fileformats/netcdf.py | 50 ++++++++++++++----- .../fileformats/netcdf/test_Saver__ugrid.py | 2 +- 2 files changed, 39 insertions(+), 13 deletions(-) diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py index 19f973191f..045f183170 100644 --- a/lib/iris/fileformats/netcdf.py +++ b/lib/iris/fileformats/netcdf.py @@ -1691,8 +1691,8 @@ def record_dimension(dim_name, length, dim_coords=[]): dimension_names.clear() mesh = cube.mesh if mesh is None: - mesh_location_dimnames = {} cube_mesh_dim = None + cube_mesh_dim_name = "" else: # Create all the mesh dimensions. # NOTE: one of these will be a cube dimension, but that is not @@ -1700,23 +1700,49 @@ def record_dimension(dim_name, length, dim_coords=[]): # (face,edge,node), and before non-mesh dimensions. mesh_location_dimnames = {} for location in MESH_LOCATIONS: - location_select_key = f"include_{location}s" - mesh_coords = mesh.coords(**{location_select_key: True}) - if len(mesh_coords) > 0: - (dim_length,) = mesh_coords[0].shape # Must be 1-d - location_dim_attr = f"{location}_dimension" - dim_name = getattr(mesh, location_dim_attr) - if dim_name is None: - dim_name = f"{mesh.name()}_{location}" - while dim_name in self._existing_dim: - self._increment_name(dim_name) + # Find if this location exists in the mesh, and a characteristic + # coordinate to identify it with. + # To use only _required_ UGRID components, we use a location + # coord for nodes, but a connectivity for faces/edges + if location == 'node': + # For definiteness, use the 'X' axis one. + dim_coords = mesh.coords(include_nodes=True, + axis='x') + else: + cf_role = f"{location}_node_connectivity" + dim_coords = mesh.connectivities(cf_role=cf_role) + if len(dim_coords) > 0: + # There should only be 1 connectivity of a given type + assert len(dim_coords) == 1 + # As the mesh contains this location, we want to include this + # dim in our returned mesh dims. + # N.B. the connectivity is used as the identifying 'dim_coord' + mesh_coord = dim_coords[0] + if mesh_coord in self._dim_coords: + # Use the name already recorded for this dim of this mesh + dim_name = self._name_coord_map.name(mesh_coord) + else: + if location == 'node': + # always 1-d + (dim_length,) = mesh_coord.shape + else: + # extract source dim, respecting dim-ordering + dim_length = mesh_coord.shape[mesh_coord.src_dim] + location_dim_attr = f"{location}_dimension" + dim_name = getattr(mesh, location_dim_attr) + if dim_name is None: + dim_name = f"{mesh.name()}_{location}" + while dim_name in self._existing_dim: + self._increment_name(dim_name) record_dimension(dim_name, dim_length) + # Store the mesh dims indexed by location mesh_location_dimnames[location] = dim_name # Identify the cube dimension which maps to the mesh. # NB always at least one mesh-coord: it has exactly one cube-dim (cube_mesh_dim,) = cube.coords(mesh_coords=True)[0].cube_dims(cube) + cube_mesh_dim_name = mesh_location_dimnames[cube.location] mesh_dimensions = dimension_names.copy() @@ -1726,7 +1752,7 @@ def record_dimension(dim_name, length, dim_coords=[]): if dim == cube_mesh_dim: # Handle a mesh dimension: we already named this. dim_coords = [] - dim_name = mesh_location_dimnames[cube.location] + dim_name = cube_mesh_dim_name else: # Get a name from the dim-coord (if any). dim_coords = cube.coords(dimensions=dim, dim_coords=True) 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 978c641269..1f34de3c19 100644 --- a/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py +++ b/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py @@ -270,7 +270,7 @@ def check_save(self, data): print(text) return scan_dataset(self.tempfile_path) - def test_basic_mesh_cdl(self): + def test_basic_mesh(self): data = make_cube(mesh=_DEFAULT_MESH) dims, vars = self.check_save(data) From 24c6697b6fd4682fce4e8d38f2f0b63aa27a248b Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 7 Sep 2021 10:34:16 +0100 Subject: [PATCH 12/26] Identify mesh-location-dims with element-coord/connectivity. --- lib/iris/fileformats/netcdf.py | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py index 045f183170..07929743e1 100644 --- a/lib/iris/fileformats/netcdf.py +++ b/lib/iris/fileformats/netcdf.py @@ -1694,10 +1694,10 @@ def record_dimension(dim_name, length, dim_coords=[]): cube_mesh_dim = None cube_mesh_dim_name = "" else: - # Create all the mesh dimensions. - # NOTE: one of these will be a cube dimension, but that is not - # special. We *do* want to list these in a priority order - # (face,edge,node), and before non-mesh dimensions. + # Identify all the mesh dimensions. + # NOTE: one of these will be a cube dimension, but that one does not + # get any special handling. We *do* want to list/create them in a + # definite order (face,edge,node), and before non-mesh dimensions. mesh_location_dimnames = {} for location in MESH_LOCATIONS: # Find if this location exists in the mesh, and a characteristic @@ -1705,21 +1705,23 @@ def record_dimension(dim_name, length, dim_coords=[]): # To use only _required_ UGRID components, we use a location # coord for nodes, but a connectivity for faces/edges if location == 'node': - # For definiteness, use the 'X' axis one. + # For nodes, identify the dim with a coordinate variable. + # Selecting the X-axis one for definiteness. dim_coords = mesh.coords(include_nodes=True, axis='x') else: + # For face/edge, use the non-optional connectivity variable. cf_role = f"{location}_node_connectivity" dim_coords = mesh.connectivities(cf_role=cf_role) if len(dim_coords) > 0: - # There should only be 1 connectivity of a given type - assert len(dim_coords) == 1 # As the mesh contains this location, we want to include this # dim in our returned mesh dims. - # N.B. the connectivity is used as the identifying 'dim_coord' + # We should have 1 identifying variable (of either type). + assert len(dim_coords) == 1 mesh_coord = dim_coords[0] if mesh_coord in self._dim_coords: - # Use the name already recorded for this dim of this mesh + # Use the previous name for this dim of this mesh, from + # a previously saved cube with the same mesh. dim_name = self._name_coord_map.name(mesh_coord) else: if location == 'node': @@ -1764,6 +1766,9 @@ def record_dimension(dim_name, length, dim_coords=[]): if coord in self._dim_coords: # Return the dim_name associated with the existing # coordinate. + # NOTE: as _name_coord_map tracks actual file variables, + # which are not created till later, this must be from + # a _previously_ saved cube. dim_name = self._name_coord_map.name(coord) else: # Determine a unique dimension name from the coord From 8f247cc820d178e6418abed50534258e6d057a8c Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 7 Sep 2021 16:59:37 +0100 Subject: [PATCH 13/26] Fix equality between connectivities of different shapes. --- lib/iris/experimental/ugrid/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/lib/iris/experimental/ugrid/__init__.py b/lib/iris/experimental/ugrid/__init__.py index bfc570fcfd..e31de8dd59 100644 --- a/lib/iris/experimental/ugrid/__init__.py +++ b/lib/iris/experimental/ugrid/__init__.py @@ -522,6 +522,8 @@ def __eq__(self, other): if hasattr(other, "metadata"): # metadata comparison eq = self.metadata == other.metadata + if eq: + eq = self.shape == other.shape if eq: eq = ( self.indices_by_src() == other.indices_by_src() From 2f4fadca9c9385aee682e8c2f87b645319b5f2c8 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 7 Sep 2021 17:00:39 +0100 Subject: [PATCH 14/26] Small test fix. --- lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py | 2 -- 1 file changed, 2 deletions(-) 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 1f34de3c19..fb0e33fad8 100644 --- a/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py +++ b/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py @@ -353,8 +353,6 @@ def test_multi_cubes_identical_meshes(self): mesh2 = make_mesh() cube1 = make_cube(var_name="a", mesh=mesh1) cube2 = make_cube(var_name="b", mesh=mesh2) - cube1 = make_cube(var_name="a") - cube2 = make_cube(var_name="b") dims, vars = self.check_save([cube1, cube2]) # there are 2 meshes in the file From 64347fd4bf3da98cafe38e60861a664433b2b7ad Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 7 Sep 2021 17:08:47 +0100 Subject: [PATCH 15/26] Tracking of mesh dims, including disambiguated: Working with existing simple tests. --- lib/iris/fileformats/netcdf.py | 61 ++++++++++++++++++---------------- 1 file changed, 32 insertions(+), 29 deletions(-) diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py index 07929743e1..9f2de92e49 100644 --- a/lib/iris/fileformats/netcdf.py +++ b/lib/iris/fileformats/netcdf.py @@ -1023,12 +1023,15 @@ def __init__(self, filename, netcdf_format): # All persistent variables #: CF name mapping with iris coordinates self._name_coord_map = CFNameCoordMap() - #: List of dimension coordinates added to the file - self._dim_coords = [] + #: Map of dimensions to characteristic coordinates with which they are identified + self._dim_names_and_coords = CFNameCoordMap() #: List of grid mappings added to the file self._coord_systems = [] #: A dictionary, listing dimension names and corresponding length self._existing_dim = {} + #: A map from meshes to their actual file dimensions (names). + # NB: might not match those of the mesh, if they were 'incremented'. + self._mesh_dims = {} #: A dictionary, mapping formula terms to owner cf variable name self._formula_terms_cache = {} #: NetCDF dataset @@ -1664,9 +1667,9 @@ def _get_dim_names(self, cube): """ dimension_names = [] - def record_dimension(dim_name, length, dim_coords=[]): + def record_dimension(dim_name, length, matching_coords=[]): """ - Record a file dimension and its length. + Record a file dimension, its length and associated coordinates. If the dimension has been seen already, check that it's length matches the earlier finding. @@ -1680,9 +1683,9 @@ def record_dimension(dim_name, length, dim_coords=[]): assert self._existing_dim[dim_name] == length # Add new coords (mesh or dim) to the already-seen list. - for coord in dim_coords: - if coord not in self._dim_coords: - self._dim_coords.append(coord) + for coord in matching_coords: + if coord not in self._dim_names_and_coords.coords: + self._dim_names_and_coords.append(dim_name, coord) # Add the latest name to the returned list dimension_names.append(dim_name) @@ -1692,7 +1695,6 @@ def record_dimension(dim_name, length, dim_coords=[]): mesh = cube.mesh if mesh is None: cube_mesh_dim = None - cube_mesh_dim_name = "" else: # Identify all the mesh dimensions. # NOTE: one of these will be a cube dimension, but that one does not @@ -1704,11 +1706,10 @@ def record_dimension(dim_name, length, dim_coords=[]): # coordinate to identify it with. # To use only _required_ UGRID components, we use a location # coord for nodes, but a connectivity for faces/edges - if location == 'node': + if location == "node": # For nodes, identify the dim with a coordinate variable. # Selecting the X-axis one for definiteness. - dim_coords = mesh.coords(include_nodes=True, - axis='x') + dim_coords = mesh.coords(include_nodes=True, axis="x") else: # For face/edge, use the non-optional connectivity variable. cf_role = f"{location}_node_connectivity" @@ -1719,12 +1720,19 @@ def record_dimension(dim_name, length, dim_coords=[]): # We should have 1 identifying variable (of either type). assert len(dim_coords) == 1 mesh_coord = dim_coords[0] - if mesh_coord in self._dim_coords: + match = mesh_coord in self._dim_names_and_coords.coords + if match: + # For mesh-identifying coords, we require the *same* + # coord, not an identical one (i.e. "is" not "==") + name = self._dim_names_and_coords.name(mesh_coord) + stored_coord = self._dim_names_and_coords.coord(name) + match = mesh_coord is stored_coord + if match: # Use the previous name for this dim of this mesh, from # a previously saved cube with the same mesh. - dim_name = self._name_coord_map.name(mesh_coord) + dim_name = self._dim_names_and_coords.name(mesh_coord) else: - if location == 'node': + if location == "node": # always 1-d (dim_length,) = mesh_coord.shape else: @@ -1734,17 +1742,18 @@ def record_dimension(dim_name, length, dim_coords=[]): dim_name = getattr(mesh, location_dim_attr) if dim_name is None: dim_name = f"{mesh.name()}_{location}" - while dim_name in self._existing_dim: - self._increment_name(dim_name) + while dim_name in self._existing_dim: + dim_name = self._increment_name(dim_name) + + record_dimension(dim_name, dim_length, dim_coords) - record_dimension(dim_name, dim_length) # Store the mesh dims indexed by location mesh_location_dimnames[location] = dim_name # Identify the cube dimension which maps to the mesh. - # NB always at least one mesh-coord: it has exactly one cube-dim (cube_mesh_dim,) = cube.coords(mesh_coords=True)[0].cube_dims(cube) - cube_mesh_dim_name = mesh_location_dimnames[cube.location] + # Record the actual file dimension names for each mesh. + self._mesh_dims[mesh] = mesh_location_dimnames.copy() mesh_dimensions = dimension_names.copy() @@ -1754,7 +1763,7 @@ def record_dimension(dim_name, length, dim_coords=[]): if dim == cube_mesh_dim: # Handle a mesh dimension: we already named this. dim_coords = [] - dim_name = cube_mesh_dim_name + 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) @@ -1763,13 +1772,10 @@ def record_dimension(dim_name, length, dim_coords=[]): coord = dim_coords[0] # always have at least one # Add only dimensions that have not already been added. - if coord in self._dim_coords: + if coord in self._dim_names_and_coords.coords: # Return the dim_name associated with the existing # coordinate. - # NOTE: as _name_coord_map tracks actual file variables, - # which are not created till later, this must be from - # a _previously_ saved cube. - dim_name = self._name_coord_map.name(coord) + dim_name = self._dim_names_and_coords.name(coord) else: # Determine a unique dimension name from the coord dim_name = self._get_coord_variable_name(cube, coord) @@ -2075,10 +2081,7 @@ def _create_mesh(self, cube, mesh): self._set_cf_var_attributes(cf_mesh_var, mesh) # Get the mesh-element dim names. - mesh_dims = { - location: getattr(mesh, f"{location}_dimension") - for location in MESH_LOCATIONS - } + mesh_dims = self._mesh_dims[mesh] # Add the element coordinate variables. for location in MESH_LOCATIONS: coords_meshobj_attr = f"{location}_coords" From 40e34906e4330f77f0a9c59ee5c90cdbd28bf02f Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Sat, 11 Sep 2021 11:08:05 +0100 Subject: [PATCH 16/26] Ugrid example #4: hybrid coord is on vertical dim *not* mesh. --- .../ugrid_conventions_examples/README.txt | 48 ++++--------------- .../ugrid_ex4_3d_layered.cdl | 8 +--- 2 files changed, 11 insertions(+), 45 deletions(-) diff --git a/lib/iris/tests/integration/experimental/ugrid_conventions_examples/README.txt b/lib/iris/tests/integration/experimental/ugrid_conventions_examples/README.txt index 2f9bcfaf5e..8221bffcb5 100644 --- a/lib/iris/tests/integration/experimental/ugrid_conventions_examples/README.txt +++ b/lib/iris/tests/integration/experimental/ugrid_conventions_examples/README.txt @@ -2,45 +2,15 @@ Examples generated from CDL example sections in UGRID conventions v1.0 ( see webpage: https://ugrid-conventions.github.io/ugrid-conventions/ ) CHANGES: - * all files had a data-var added, for ease of iris-roundtripping + * added a data-var to all examples, for ease of iris-roundtripping * EX4 : - had a couple of missing ";"s at lineends - - existing var "Mesh2_surface" is tied to the node dimension but has ':location = "face"' - - actually, ALL the formula terms should map to 'face', not nodes. - - created data-var maps (layers, faces), as reqd - -/home/h05/itpp/git/iris/iris_main/lib/iris/tests/results/ugrid_ref - $ for f in $(ls *.cdl); do n=$(echo $f | grep -o "[^.]*" | grep "ugrid"); echo $n; ncgen $n.cdl -4 -o $n.nc; done - ugrid_ex1_1d_mesh - ugrid_ex2_2d_triangular - ugrid_ex3_2d_flexible - ugrid_ex4_3d_layered - - (ncdump -h $n.nc >$n.ncdump.txt) - $ for f in $(ls *.cdl); do n=$(echo $f | grep -o "[^.]*" | grep "ugrid"); echo $n; ncdump -h $n.nc >$n.ncdump.txt; done - - (xxdiff $n.cdl $n.ncdump.txt;) - $ for f in $(ls *.cdl); do n=$(echo $f | grep -o "[^.]*" | grep "ugrid"); echo $n; xxdiff $n.cdl $n.ncdump.txt; done - -Then for compatibility testing... - (python iris_loadsave.py $n.nc) - $ for f in $(ls *.cdl); do n=$(echo $f | grep -o "[^.]*" | grep "ugrid"); echo $n; python iris_loadsave.py $n.nc; done - - - -So each of (4) examples has : - .cdl : original text from UGRID webpage - .nc : "ncgen" output = generated netcdf file - .ncdump.txt : "ncdump -h" output = re-generated CDL from nc file -(from iris_loadsave.py) - _REDUMP_cdl.txt : same as .ncdump.txt - _RESAVED.nc : from loading .nc + re-saving it - _RESAVED_REDUMP_cdl.txt : ncdump of _RESAVED.nc - -NEWSTYLE CHECKS ? -.CDL : - ==> ncgen ==> .nc - ==> load ==> ex_iris_data - ==> save ==> ex_iris_resaved - ==> load ==> ex_iris_saveload :: COMPARE with ex_iris_data + - the formula terms (depth+surface) should map to 'Mesh2_layers', and not to the mesh at all. + - use Mesh2d_layers dim, and have no 'mesh' or 'location' + * "EX4a" -- closer mix of hybrid-vertical and mesh dimensions + - *don't* think we can have a hybrid coord ON the mesh dimension + - mesh being a vertical location (only) seems to make no sense + - .. and implies that the mesh is 1d and ordered, which is not really unstructured at all + - *could* have hybrid-height with the _orography_ mapping to the mesh + - doesn't match the UGRID examples, but see : iris.tests.unit.fileformats.netcdf.test_Saver__ugrid.TestSaveUgrid__cube.test_nonmesh_hybrid_dim diff --git a/lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex4_3d_layered.cdl b/lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex4_3d_layered.cdl index 5b3592c5d8..d154502018 100644 --- a/lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex4_3d_layered.cdl +++ b/lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex4_3d_layered.cdl @@ -100,18 +100,14 @@ Mesh2_layers:standard_name = "ocean_sigma_coordinate" ; Mesh2_layers:long_name = "sigma at layer midpoints" ; Mesh2_layers:positive = "up" ; Mesh2_layers:formula_terms = "sigma: Mesh2_layers eta: Mesh2_surface depth: Mesh2_depth" ; -double Mesh2_depth(nMesh2_face) ; +double Mesh2_depth(Mesh2_layers) ; Mesh2_depth:standard_name = "sea_floor_depth_below_geoid" ; Mesh2_depth:units = "m" ; Mesh2_depth:positive = "down" ; -Mesh2_depth:mesh = "Mesh2" ; -Mesh2_depth:location = "face" ; Mesh2_depth:coordinates = "Mesh2_node_x Mesh2_node_y" ; -double Mesh2_surface(nMesh2_face) ; +double Mesh2_surface(Mesh2_layers) ; Mesh2_surface:standard_name = "sea_surface_height_above_geoid" ; Mesh2_surface:units = "m" ; -Mesh2_surface:mesh = "Mesh2" ; -Mesh2_surface:location = "face" ; Mesh2_surface:coordinates = "Mesh2_face_x Mesh2_face_y" ; float datavar(Mesh2_layers, nMesh2_face) ; From 1d5e3d9f72c70a1d9a8349a3853341d7f1b93bc5 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Mon, 13 Sep 2021 09:41:19 +0100 Subject: [PATCH 17/26] Avoid non-ugrid problem in integration roundtrip 'ex4' test. --- .../experimental/test_ugrid_save.py | 11 ++++++++--- .../TestBasicSave/ugrid_ex4_3d_layered.cdl | 18 ++++++------------ 2 files changed, 14 insertions(+), 15 deletions(-) diff --git a/lib/iris/tests/integration/experimental/test_ugrid_save.py b/lib/iris/tests/integration/experimental/test_ugrid_save.py index 236ba006b0..13ef3956fd 100644 --- a/lib/iris/tests/integration/experimental/test_ugrid_save.py +++ b/lib/iris/tests/integration/experimental/test_ugrid_save.py @@ -73,9 +73,6 @@ def test_example_roundtrips(self): # for data derived from each UGRID example CDL. for ex_name, filepath in self.example_names_paths.items(): print(f"Roundtrip checking : {ex_name}") - if "ex4" in ex_name: - # Skip this one now : still causing some problems ... - continue target_ncfile_path = str(self.temp_dir / f"{ex_name}.nc") # Create a netcdf file from the test CDL. check_call( @@ -86,11 +83,18 @@ def test_example_roundtrips(self): # Load the original as Iris data with PARSE_UGRID_ON_LOAD.context(): orig_cubes = iris.load(target_ncfile_path) + + if "ex4" in ex_name: + # Discard the extra formula terms component cubes + # Saving these does not do what you expect + orig_cubes = orig_cubes.extract('datavar') + # Save-and-load-back to compare the Iris saved result. resave_ncfile_path = str(self.temp_dir / f"{ex_name}_resaved.nc") iris.save(orig_cubes, resave_ncfile_path) with PARSE_UGRID_ON_LOAD.context(): savedloaded_cubes = iris.load(resave_ncfile_path) + # This should match the original exactly # ..EXCEPT for our inability to compare meshes. for orig, reloaded in zip(orig_cubes, savedloaded_cubes): @@ -99,6 +103,7 @@ def test_example_roundtrips(self): cube.attributes.pop("Conventions", None) # Remove var-names, which may differ. cube.var_name = None + # Compare the mesh contents (as we can't compare actual meshes) self.assertEqual(orig.location, reloaded.location) orig_mesh = orig.mesh diff --git a/lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex4_3d_layered.cdl b/lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex4_3d_layered.cdl index 05517f340f..64962e79aa 100644 --- a/lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex4_3d_layered.cdl +++ b/lib/iris/tests/results/integration/experimental/ugrid_save/TestBasicSave/ugrid_ex4_3d_layered.cdl @@ -78,28 +78,22 @@ variables: Mesh2_layers:long_name = "sigma at layer midpoints" ; Mesh2_layers:positive = "up" ; Mesh2_layers:formula_terms = "sigma: Mesh2_layers eta: Mesh2_surface depth: Mesh2_depth" ; - double Mesh2_depth(nMesh2_face) ; + double Mesh2_depth(Mesh2_layers) ; Mesh2_depth:units = "m" ; Mesh2_depth:standard_name = "sea_floor_depth_below_geoid" ; - Mesh2_depth:location = "face" ; - Mesh2_depth:mesh = "Mesh2" ; Mesh2_depth:positive = "down" ; - double Mesh2_surface(nMesh2_face) ; + double Mesh2_surface(Mesh2_layers) ; Mesh2_surface:units = "m" ; Mesh2_surface:standard_name = "sea_surface_height_above_geoid" ; - Mesh2_surface:location = "face" ; - Mesh2_surface:mesh = "Mesh2" ; - double Mesh2_depth_0(nMesh2_face) ; + double Mesh2_depth_0(Mesh2_layers) ; Mesh2_depth_0:standard_name = "sea_floor_depth_below_geoid" ; Mesh2_depth_0:units = "m" ; Mesh2_depth_0:positive = "down" ; - Mesh2_depth_0:mesh = "Mesh2" ; - Mesh2_depth_0:location = "face" ; - double Mesh2_surface_0(nMesh2_face) ; + Mesh2_depth_0:coordinates = "Mesh2_depth Mesh2_surface" ; + double Mesh2_surface_0(Mesh2_layers) ; Mesh2_surface_0:standard_name = "sea_surface_height_above_geoid" ; Mesh2_surface_0:units = "m" ; - Mesh2_surface_0:mesh = "Mesh2" ; - Mesh2_surface_0:location = "face" ; + Mesh2_surface_0:coordinates = "Mesh2_depth Mesh2_surface" ; // global attributes: :Conventions = "CF-1.7" ; From 1c0d36117d4500b44050da00b17c75de8149c7b8 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 14 Sep 2021 00:50:36 +0100 Subject: [PATCH 18/26] Fix example files readme. --- .../experimental/ugrid_conventions_examples/README.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/iris/tests/integration/experimental/ugrid_conventions_examples/README.txt b/lib/iris/tests/integration/experimental/ugrid_conventions_examples/README.txt index 8221bffcb5..2a9b5bde35 100644 --- a/lib/iris/tests/integration/experimental/ugrid_conventions_examples/README.txt +++ b/lib/iris/tests/integration/experimental/ugrid_conventions_examples/README.txt @@ -7,7 +7,7 @@ CHANGES: - had a couple of missing ";"s at lineends - the formula terms (depth+surface) should map to 'Mesh2_layers', and not to the mesh at all. - use Mesh2d_layers dim, and have no 'mesh' or 'location' - * "EX4a" -- closer mix of hybrid-vertical and mesh dimensions + * "EX4a" -- possibly (future) closer mix of hybrid-vertical and mesh dimensions - *don't* think we can have a hybrid coord ON the mesh dimension - mesh being a vertical location (only) seems to make no sense - .. and implies that the mesh is 1d and ordered, which is not really unstructured at all From c462a8369386dae8bc73381c36d2fc6f79fa63cb Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 14 Sep 2021 00:51:51 +0100 Subject: [PATCH 19/26] More tests - working. --- lib/iris/fileformats/netcdf.py | 20 +- .../TestSaveUgrid__cube/basic_mesh.cdl | 29 ++ .../fileformats/netcdf/test_Saver__ugrid.py | 276 +++++++++++++++--- 3 files changed, 284 insertions(+), 41 deletions(-) create mode 100644 lib/iris/tests/results/unit/fileformats/netcdf/Saver__ugrid/TestSaveUgrid__cube/basic_mesh.cdl diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py index 9f2de92e49..fe636d1a49 100644 --- a/lib/iris/fileformats/netcdf.py +++ b/lib/iris/fileformats/netcdf.py @@ -2114,21 +2114,27 @@ def _create_mesh(self, cube, mesh): if last_dim not in self._dataset.dimensions: length = conn.shape[1] self._dataset.createDimension(last_dim, length) + # Create variable. + # NOTE: for connectivities *with missing points*, this will use a + # fixed standard fill-value of -1, and add a _FillValue property. conn_dims = (mesh_dims[loc_from], last_dim) cf_conn_name = self._create_generic_cf_array_var( - cube, [], conn, element_dims=conn_dims + cube, [], conn, element_dims=conn_dims, fill_value=-1 ) # Add essential attributes to the Connectivity variable. cf_conn_var = self._dataset.variables[cf_conn_name] _setncattr(cf_conn_var, "cf_role", cf_conn_attr_name) _setncattr(cf_conn_var, "start_index", conn.start_index) + # If content was masked, also add a "_FillValue" property. + # N.B. for now, we are using -1 as a universal 'safe' value. + if np.ma.is_masked(conn.indices): + _setncattr(cf_var, "_FillValue", -1) + # Record the connectivity on the parent mesh var. _setncattr(cf_mesh_var, cf_conn_attr_name, cf_conn_name) - # TODO: possibly we need to handle boundary coords/conns separately - return cf_mesh_name def _set_cf_var_attributes(self, cf_var, element): @@ -2170,7 +2176,7 @@ 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 + self, cube, cube_dim_names, element, element_dims=None, fill_value=None ): """ Create the associated CF-netCDF variable in the netCDF dataset for the @@ -2196,6 +2202,8 @@ def _create_generic_cf_array_var( otherwise these are taken from `element.cube_dims[cube]`. For Mesh components (element coordinates and connectivities), this *must* be passed in, as "element.cube_dims" does not function. + * fill_value (number or None): + If set, fill any masked data points with this value. Returns: var_name (string): @@ -2258,6 +2266,10 @@ def _create_generic_cf_array_var( element_type = type(element).__name__ data = self._ensure_valid_dtype(data, element_type, element) + if fill_value is not None: + # Use a specific fill-value in place of the netcdf default. + data = np.ma.filled(data, fill_value) + # Check if this is a dim-coord. is_dimcoord = element in cube.dim_coords diff --git a/lib/iris/tests/results/unit/fileformats/netcdf/Saver__ugrid/TestSaveUgrid__cube/basic_mesh.cdl b/lib/iris/tests/results/unit/fileformats/netcdf/Saver__ugrid/TestSaveUgrid__cube/basic_mesh.cdl new file mode 100644 index 0000000000..91516ddae3 --- /dev/null +++ b/lib/iris/tests/results/unit/fileformats/netcdf/Saver__ugrid/TestSaveUgrid__cube/basic_mesh.cdl @@ -0,0 +1,29 @@ +dimensions: + Mesh2d_face_N_nodes = 4 ; + Mesh2d_faces = 2 ; + Mesh2d_nodes = 5 ; +variables: + int Mesh2d ; + Mesh2d:cf_role = "mesh_topology" ; + Mesh2d:topology_dimension = 2 ; + Mesh2d:node_coordinates = "node_x node_y" ; + Mesh2d:face_coordinates = "face_x face_y" ; + Mesh2d:face_node_connectivity = "mesh2d_faces" ; + int64 node_x(Mesh2d_nodes) ; + node_x:standard_name = "longitude" ; + int64 node_y(Mesh2d_nodes) ; + node_y:standard_name = "latitude" ; + int64 face_x(Mesh2d_faces) ; + face_x:standard_name = "longitude" ; + int64 face_y(Mesh2d_faces) ; + face_y:standard_name = "latitude" ; + int mesh2d_faces(Mesh2d_faces, Mesh2d_face_N_nodes) ; + mesh2d_faces:cf_role = "face_node_connectivity" ; + mesh2d_faces:start_index = 0LL ; + float unknown(Mesh2d_faces) ; + unknown:mesh = "Mesh2d" ; + unknown:location = "face" ; + +// global attributes: + :Conventions = "CF-1.7" ; +} 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 fb0e33fad8..7dcd8b0f9c 100644 --- a/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py +++ b/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py @@ -14,12 +14,16 @@ from subprocess import check_output import tempfile +import netCDF4 as nc import numpy as np from iris import save from iris.coords import AuxCoord -from iris.cube import Cube +from iris.cube import Cube, CubeList from iris.experimental.ugrid import Connectivity, Mesh +from iris.tests.stock import realistic_4d +from iris.util import new_axis, demote_dim_coord_to_aux_coord + XY_LOCS = ("x", "y") XY_NAMES = ("longitude", "latitude") @@ -110,6 +114,8 @@ def apply_xyargs(coords, xyargs): connectivities=connectivities.values(), ) applyargs(mesh, mesh_kwargs) + + # N.B. *nasty!* attach an extra convenience property for handling purposes. mesh._mesh_dims = mesh_dims return mesh @@ -155,6 +161,19 @@ def make_cube(mesh=_DEFAULT_MESH, location="face", **kwargs): return cube +def add_height_dim(cube): + # Add an extra inital dimension to the cube. + cube = cube.copy() # Avoid trashing the input cube. + cube.add_aux_coord(AuxCoord([0.], standard_name='height', units='m')) + # Make three copies with different heights + cubes = [cube.copy() for _ in range(3)] + for i_cube, cube in enumerate(cubes): + cube.coord('height').points = [i_cube] + # Merge to create an additional 'height' dimension. + cube = CubeList(cubes).merge_cube() + return cube + + def scan_dataset(filepath): """ Snapshot a netcdf dataset (the key metadata). @@ -166,12 +185,10 @@ def scan_dataset(filepath): * varsdict (dict): A map of each variable's properties, {var_name: propsdict} Each propsdict is {attribute-name: value} over the var's ncattrs(). - Each propsdict ALSO contains a ['_dims'] entry listing the + Each propsdict ALSO contains a ['_DIMS'] entry listing the variable's dims. """ - import netCDF4 as nc - ds = nc.Dataset(filepath) # dims dict is {name: len} dimsdict = {name: dim.size for name, dim in ds.dimensions.items()} @@ -179,7 +196,7 @@ def scan_dataset(filepath): varsdict = {} for name, var in ds.variables.items(): varsdict[name] = {prop: getattr(var, prop) for prop in var.ncattrs()} - varsdict[name]["_dims"] = list(var.dimensions) + varsdict[name]["_DIMS"] = list(var.dimensions) ds.close() return dimsdict, varsdict @@ -216,14 +233,14 @@ def vars_w_dims(varsdict, dim_names): varsdict = { name: propsdict for name, propsdict in varsdict.items() - if all(dim in propsdict["_dims"] for dim in dim_names) + if all(dim in propsdict["_DIMS"] for dim in dim_names) } return varsdict -def vars_meshvars(vars): - """Subset a varsdict, returning those which are mesh variables (by cf_role).""" - return vars_w_props(vars, cf_role="mesh_topology") +def vars_meshnames(vars): + """Return the names of all the mesh variables (found by cf_role).""" + return list(vars_w_props(vars, cf_role="mesh_topology").keys()) def vars_meshdim(vars, location, mesh_name=None): @@ -243,13 +260,16 @@ def vars_meshdim(vars, location, mesh_name=None): dim_name (string) The dim-name of the mesh dim for the given location. + TODO: relies on the element having coordinates, which in future will not + always be the case. This can be fixed + """ if mesh_name is None: # Find "the" meshvar -- assuming there is just one. - (mesh_name,) = vars_meshvars(vars).keys() + (mesh_name,) = vars_meshnames(vars) mesh_props = vars[mesh_name] loc_coords = mesh_props[f"{location}_coordinates"].split(" ") - (single_location_dim,) = vars[loc_coords[0]]["_dims"] + (single_location_dim,) = vars[loc_coords[0]]["_DIMS"] return single_location_dim @@ -271,15 +291,19 @@ def check_save(self, data): return scan_dataset(self.tempfile_path) def test_basic_mesh(self): - data = make_cube(mesh=_DEFAULT_MESH) + # Save a small mesh example and check most aspects of the resultin file. + data = make_cube() # A simple face-mapped data example. + + # Save and snapshot the result dims, vars = self.check_save(data) # There is exactly 1 mesh var. - ((mesh_name, mesh_props),) = vars_meshvars(vars).items() + (mesh_name,) = vars_meshnames(vars) # There is exactly 1 mesh-linked (data)var data_vars = vars_w_props(vars, mesh="*") ((a_name, a_props),) = data_vars.items() + mesh_props = vars[mesh_name] # The mesh var links to the mesh, with location 'faces' self.assertEqual(a_name, "unknown") @@ -289,51 +313,77 @@ def test_basic_mesh(self): # There are 2 face coords == those listed in the mesh face_coords = mesh_props["face_coordinates"].split(" ") self.assertEqual(len(face_coords), 2) - # get face dim (actually, from the dims of the first face coord) - face_dim = vars_meshdim(vars, "face") + # The face coords should both map that single dim. + face_dim = vars_meshdim(vars, "face") self.assertTrue( - all(vars[co]["_dims"] == [face_dim] for co in face_coords) + all(vars[co]["_DIMS"] == [face_dim] for co in face_coords) ) # The dims of the datavar also == [] - self.assertEqual(a_props["_dims"], [face_dim]) + self.assertEqual(a_props["_DIMS"], [face_dim]) # There are 2 node coordinates == those listed in the mesh. node_coords = mesh_props["node_coordinates"].split(" ") self.assertEqual(len(node_coords), 2) - # These should be the only ones using the 'nodes' dimension. + # These are the *only* ones using the 'nodes' dimension. node_dim = vars_meshdim(vars, "node") self.assertEqual( sorted(node_coords), sorted(vars_w_dims(vars, [node_dim]).keys()) ) # There are no edges. + self.assertNotIn('edge_node_connectivity', mesh_props) self.assertEqual( len(vars_w_props(vars, cf_role="edge_node_connectivity")), 0 ) + # The dims are precisely (nodes, faces, nodes-per-face), in that order. + self.assertEqual( + list(dims.keys()), + ['Mesh2d_nodes', 'Mesh2d_faces', 'Mesh2d_face_N_nodes']) + + # The variables are (mesh, 2*node-coords, 2*face-coords, face-nodes, data), + # in that order + self.assertEqual( + list(vars.keys()), + ['Mesh2d', + 'node_x', 'node_y', + 'face_x', 'face_y', + 'mesh2d_faces', + 'unknown', + ] + ) + + # For definiteness, also check against a full CDL snapshot + self.assertCDL(self.tempfile_path) + def test_multi_cubes_common_mesh(self): cube1 = make_cube(var_name="a") cube2 = make_cube(var_name="b") + + # Save and snapshot the result dims, vars = self.check_save([cube1, cube2]) + # there is exactly 1 mesh in the file - ((mesh_name, _),) = vars_meshvars(vars).items() + (mesh_name,) = vars_meshnames(vars) + # both the main variables reference the same mesh, and 'face' location v_a, v_b = vars["a"], vars["b"] self.assertEqual(v_a["mesh"], mesh_name) self.assertEqual(v_a["location"], "face") self.assertEqual(v_b["mesh"], mesh_name) self.assertEqual(v_b["location"], "face") - # self.assertCDL(self.tempfile_path) def test_multi_cubes_different_locations(self): cube1 = make_cube(var_name="a", location="face") cube2 = make_cube(var_name="b", location="node") + + # Save and snapshot the result dims, vars = self.check_save([cube1, cube2]) # there is exactly 1 mesh in the file - ((mesh_name, mesh_props),) = vars_meshvars(vars).items() + (mesh_name,) = vars_meshnames(vars) # the main variables reference the same mesh at different locations v_a, v_b = vars["a"], vars["b"] @@ -345,44 +395,196 @@ def test_multi_cubes_different_locations(self): # the main variables map the face and node dimensions face_dim = vars_meshdim(vars, "face") node_dim = vars_meshdim(vars, "node") - self.assertEqual(v_a["_dims"], [face_dim]) - self.assertEqual(v_b["_dims"], [node_dim]) + self.assertEqual(v_a["_DIMS"], [face_dim]) + self.assertEqual(v_b["_DIMS"], [node_dim]) def test_multi_cubes_identical_meshes(self): + # Make 2 identical meshes + # NOTE: *can't* name these explicitly, as it stops them being identical. mesh1 = make_mesh() mesh2 = make_mesh() cube1 = make_cube(var_name="a", mesh=mesh1) cube2 = make_cube(var_name="b", mesh=mesh2) - dims, vars = self.check_save([cube1, cube2]) - # there are 2 meshes in the file - mesh_vars = vars_meshvars(vars) - self.assertEqual(len(mesh_vars), 2) + # Save and snapshot the result + dims, vars = self.check_save([cube1, cube2]) - # there are two (data)variables with a 'mesh' property + # 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 + self.assertEqual(vars_meshdim(vars, 'node', mesh_name='Mesh2d'), + 'Mesh2d_nodes') + self.assertEqual(vars_meshdim(vars, 'face', mesh_name='Mesh2d'), + 'Mesh2d_faces') + 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') + + # there are exactly two data-variables with a 'mesh' property mesh_datavars = vars_w_props(vars, mesh="*") - self.assertEqual(2, len(mesh_datavars)) - self.assertEqual(["a", "b"], sorted(mesh_datavars.keys())) + self.assertEqual(["a", "b"], list(mesh_datavars)) # the data variables reference the two separate meshes a_props, b_props = vars["a"], vars["b"] - mesh_a, mesh_b = a_props["mesh"], b_props["mesh"] - self.assertEqual(sorted([mesh_a, mesh_b]), sorted(mesh_vars.keys())) + self.assertEqual(a_props['mesh'], 'Mesh2d') + self.assertEqual(a_props['location'], 'face') + self.assertEqual(b_props['mesh'], 'Mesh2d_0') + self.assertEqual(b_props['location'], 'face') + + # the data variables map the appropriate node dimensions + self.assertEqual(a_props["_DIMS"], ['Mesh2d_faces']) + self.assertEqual(b_props["_DIMS"], ['Mesh2d_faces_0']) + def test_multi_cubes_different_mesh(self): + # Check that we can correctly distinguish 2 different meshes. cube1 = make_cube(var_name="a") cube2 = make_cube(var_name="b", mesh=make_mesh(n_faces=4)) - self.check_save([cube1, cube2]) + + # Save and snapshot the result dims, vars = self.check_save([cube1, cube2]) + # there are 2 meshes in the file - mesh_vars = vars_w_props(vars, cf_role="mesh_topology") - self.assertEqual(len(mesh_vars), 2) + mesh_names = vars_meshnames(vars) + self.assertEqual(len(mesh_names), 2) + # there are two (data)variables with a 'mesh' property mesh_datavars = vars_w_props(vars, mesh="*") self.assertEqual(2, len(mesh_datavars)) self.assertEqual(["a", "b"], sorted(mesh_datavars.keys())) - # the main variables reference the same mesh, and 'face' location - # self.assertCDL(self.tempfile_path) + + # the main variables reference the respective meshes, and 'face' location + a_props, b_props = vars["a"], vars["b"] + mesh_a, loc_a = a_props['mesh'], a_props['location'] + mesh_b, loc_b = b_props['mesh'], b_props['location'] + self.assertNotEqual(mesh_a, mesh_b) + self.assertEqual(loc_a, 'face') + self.assertEqual(loc_b, 'face') + + def test_nonmesh_dim(self): + # Check where the data variable has a 'normal' dim and a mesh dim. + cube = make_cube() + cube = add_height_dim(cube) + + # Save and snapshot the result + dims, vars = self.check_save(cube) + + # have just 1 mesh, including a face and node coordinates. + (mesh_name,) = vars_meshnames(vars) + face_dim = vars_meshdim(vars, 'face', mesh_name) + _ = vars_meshdim(vars, 'node', mesh_name) + + # have just 1 data-variable + ((data_name, data_props),) = vars_w_props(vars, mesh='*').items() + + # data maps to the height + mesh dims + self.assertEqual(data_props['_DIMS'], ['height', face_dim]) + self.assertEqual(data_props['mesh'], mesh_name) + self.assertEqual(data_props['location'], 'face') + + def test_nonmesh_hybrid_dim(self): + # Check a case with a hybrid non-mesh dimension + cube = realistic_4d() + # Strip off the time and longtude dims, to make it simpler. + cube = cube[0, ..., 0] + # Remove all the unwanted coords (also loses the coord-system) + lose_coords = ('time', 'forecast_period', + 'grid_longitude', 'grid_latitude') + for coord in lose_coords: + cube.remove_coord(coord) + + # Add a mesh on the remaining (now anonymous) horizontal dimension. + i_horizontal_dim = len(cube.shape) - 1 + n_places = cube.shape[i_horizontal_dim] + mesh = make_mesh( + n_faces=n_places, + n_nodes=30 # arbitrary + unrealistic, but doesn't actually matter + ) + # Attach the mesh by adding MeshCoords + for coord in mesh.to_MeshCoords('face'): + cube.add_aux_coord(coord, (i_horizontal_dim,)) + + # Save and snapshot the result + dims, vars = self.check_save(cube) + + # have just 1 mesh, including face and node coordinates. + (mesh_name,) = vars_meshnames(vars) + face_dim = vars_meshdim(vars, 'face', mesh_name) + _ = vars_meshdim(vars, 'node', mesh_name) + + # have hybrid vertical dimension, with all the usual term variables. + self.assertIn('model_level_number', dims) + vert_vars = list(vars_w_dims(vars, ['model_level_number']).keys()) + # The list of file variables mapping the vertical dimensio: + # = the data-var, plus all the height terms + self.assertEqual( + vert_vars, + ['air_potential_temperature', 'model_level_number', + 'level_height', 'level_height_bnds', 'sigma', 'sigma_bnds'] + ) + + # have just 1 data-variable, which maps to hybrid-height and mesh dims + ((data_name, data_props),) = vars_w_props(vars, mesh='*').items() + self.assertEqual(data_props['_DIMS'], ['model_level_number', face_dim]) + self.assertEqual(data_props['mesh'], mesh_name) + self.assertEqual(data_props['location'], 'face') + + def test_alternate_dim_order(self): + # A cube transposed from the 'usual' order + # Should work much the same as the "basic" case. + cube_1 = make_cube(var_name='a') + cube_1 = add_height_dim(cube_1) + + cube_2 = cube_1.copy() + cube_2.var_name = 'b' + cube_2.transpose() + + # Save and snapshot the result + dims, vars = self.check_save([cube_1, cube_2]) + + # There is only 1 mesh + (mesh_name,) = vars_meshnames(vars) + + # both variables reference the same mesh + v_a, v_b = vars['a'], vars['b'] + self.assertEqual(v_a['mesh'], mesh_name) + self.assertEqual(v_a['location'], 'face') + self.assertEqual(v_b['mesh'], mesh_name) + self.assertEqual(v_b['location'], 'face') + + # Check the var dimensions + self.assertEqual(v_a['_DIMS'], ['height', 'Mesh2d_faces']) + self.assertEqual(v_b['_DIMS'], ['Mesh2d_faces', 'height']) + + # def test_mesh_nonstandard_dims(self): + # # A mesh with connectivities in the 'wrong' order + # # This should record the relevant dimension + # mesh = make_mesh() + # # Get the face-node connectivity + # face_nodes_conn = mesh.face_node_connectivity + # # Transpose it : N.B. this sets src_dim=1, asit should be. + # facenodes_nodeface_conn = face_nodes_conn.transpose() + # # Make a new mesh, based on this. + # mesh2 = Mesh( + # topology_dimension=mesh.topology_dimension, + # node_coords_and_axes=zip(mesh.node_coords, XY_LOCS), + # face_coords_and_axes=zip(mesh.face_coords, XY_LOCS), + # connectivities=[facenodes_nodeface_conn], + # ) + # # N.B. replicate our nasty extra testing property (so we can make_cube) + # mesh2._mesh_dims = mesh._mesh_dims + # + # # Build a cube on the modified mesh + # cube = make_cube(mesh=mesh2) + # + # # Save and snapshot the result + # dims, vars = self.check_save(cube) + # + # t_dbg = 0 + # pass if __name__ == "__main__": From 340187b07416c33e2cb38d22920bf8063176b8db Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 14 Sep 2021 00:52:36 +0100 Subject: [PATCH 20/26] Interim fix to cube summary. --- lib/iris/_representation/cube_summary.py | 11 ++++++++++- .../cube_printout/test_CubePrintout.py | 18 ++++++++++++++++++ .../cube_summary/test_CubeSummary.py | 15 ++++++++++++++- 3 files changed, 42 insertions(+), 2 deletions(-) diff --git a/lib/iris/_representation/cube_summary.py b/lib/iris/_representation/cube_summary.py index c7d0e15e59..3f049c879c 100644 --- a/lib/iris/_representation/cube_summary.py +++ b/lib/iris/_representation/cube_summary.py @@ -259,8 +259,16 @@ def __init__(self, cube, shorten=False, name_padding=35): vector_dim_coords = [ coord for coord in dim_coords if id(coord) not in scalar_coord_ids ] + if cube.mesh is None: + mesh_coords = [] + else: + mesh_coords = [coord for coord in aux_coords + if hasattr(coord, 'mesh')] + vector_aux_coords = [ - coord for coord in aux_coords if id(coord) not in scalar_coord_ids + coord for coord in aux_coords + if (id(coord) not in scalar_coord_ids + and coord not in mesh_coords) ] vector_derived_coords = [ coord @@ -300,6 +308,7 @@ def add_vector_section(title, contents, iscoord=True): ) add_vector_section("Dimension coordinates:", vector_dim_coords) + add_vector_section("Mesh coordinates:", mesh_coords) add_vector_section("Auxiliary coordinates:", vector_aux_coords) add_vector_section("Derived coordinates:", vector_derived_coords) add_vector_section("Cell measures:", vector_cell_measures, False) diff --git a/lib/iris/tests/unit/representation/cube_printout/test_CubePrintout.py b/lib/iris/tests/unit/representation/cube_printout/test_CubePrintout.py index f49c9f9c0c..8370c719f0 100644 --- a/lib/iris/tests/unit/representation/cube_printout/test_CubePrintout.py +++ b/lib/iris/tests/unit/representation/cube_printout/test_CubePrintout.py @@ -18,6 +18,7 @@ DimCoord, ) from iris.cube import Cube +from iris.tests.stock.mesh import sample_mesh_cube class TestCubePrintout___str__(tests.IrisTest): @@ -513,6 +514,23 @@ def test_section_cell_methods(self): ] self.assertEqual(rep, expected) + def test_unstructured_cube(self): + # Check a sample mesh-cube against the expected result. + cube = sample_mesh_cube() + rep = cube_replines(cube) + expected = [ + "mesh_phenom / (unknown) (level: 2; i_mesh_face: 3)", + " Dimension coordinates:", + " level x -", + " i_mesh_face - x", + " Mesh coordinates:", + " latitude - x", + " longitude - x", + " Auxiliary coordinates:", + " mesh_face_aux - x", + ] + self.assertEqual(rep, expected) + if __name__ == "__main__": tests.main() diff --git a/lib/iris/tests/unit/representation/cube_summary/test_CubeSummary.py b/lib/iris/tests/unit/representation/cube_summary/test_CubeSummary.py index 79baf65c8b..c8af3437e6 100644 --- a/lib/iris/tests/unit/representation/cube_summary/test_CubeSummary.py +++ b/lib/iris/tests/unit/representation/cube_summary/test_CubeSummary.py @@ -20,6 +20,7 @@ DimCoord, ) from iris.cube import Cube +from iris.tests.stock.mesh import sample_mesh_cube def example_cube(): @@ -56,6 +57,7 @@ def test_blank_cube(self): expected_vector_sections = [ "Dimension coordinates:", + "Mesh coordinates:", "Auxiliary coordinates:", "Derived coordinates:", "Cell measures:", @@ -211,7 +213,7 @@ def test_scalar_cube(self): self.assertEqual(rep.header.dimension_header.dim_names, []) self.assertEqual(rep.header.dimension_header.shape, []) self.assertEqual(rep.header.dimension_header.contents, ["scalar cube"]) - self.assertEqual(len(rep.vector_sections), 5) + self.assertEqual(len(rep.vector_sections), 6) self.assertTrue( all(sect.is_empty() for sect in rep.vector_sections.values()) ) @@ -303,6 +305,17 @@ def test_attributes_subtle_differences(self): self.assertEqual(co3a_summ.extra, "arr1=array([5, 6])") self.assertEqual(co3b_summ.extra, "arr1=array([[5], [6]])") + def test_unstructured_cube(self): + cube = sample_mesh_cube() + rep = CubeSummary(cube) + # Just check that coordinates appear in the expected sections + dim_section = rep.vector_sections["Dimension coordinates:"] + mesh_section = rep.vector_sections["Mesh coordinates:"] + aux_section = rep.vector_sections["Auxiliary coordinates:"] + self.assertEqual(len(dim_section.contents), 2) + self.assertEqual(len(mesh_section.contents), 2) + self.assertEqual(len(aux_section.contents), 1) + if __name__ == "__main__": tests.main() From c44f8fe933a24227595717fae1b34fab9b7b2f0d Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 14 Sep 2021 00:52:58 +0100 Subject: [PATCH 21/26] Interim fix connectivity comparison. --- lib/iris/experimental/ugrid/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/lib/iris/experimental/ugrid/__init__.py b/lib/iris/experimental/ugrid/__init__.py index e31de8dd59..91e5359cc0 100644 --- a/lib/iris/experimental/ugrid/__init__.py +++ b/lib/iris/experimental/ugrid/__init__.py @@ -3977,6 +3977,8 @@ def _build_mesh_coords(mesh, cf_var): } mesh_dim_name = locations_dimensions[cf_var.location] # (Only expecting 1 mesh dimension per cf_var). + if mesh_dim_name not in cf_var.dimensions: + oops = True mesh_dim = cf_var.dimensions.index(mesh_dim_name) mesh_coords = mesh.to_MeshCoords(location=cf_var.location) From 7e0b3765fa2a791920f97619438904ad67ba288d Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 14 Sep 2021 11:28:17 +0100 Subject: [PATCH 22/26] Support connectivity missing indices and alternate dims - no test for missing yet. --- lib/iris/_representation/cube_summary.py | 11 +- lib/iris/fileformats/netcdf.py | 22 +- .../experimental/test_ugrid_save.py | 2 +- .../fileformats/netcdf/test_Saver__ugrid.py | 240 +++++++++++------- 4 files changed, 169 insertions(+), 106 deletions(-) diff --git a/lib/iris/_representation/cube_summary.py b/lib/iris/_representation/cube_summary.py index 3f049c879c..68f86832f5 100644 --- a/lib/iris/_representation/cube_summary.py +++ b/lib/iris/_representation/cube_summary.py @@ -262,13 +262,14 @@ def __init__(self, cube, shorten=False, name_padding=35): if cube.mesh is None: mesh_coords = [] else: - mesh_coords = [coord for coord in aux_coords - if hasattr(coord, 'mesh')] + mesh_coords = [ + coord for coord in aux_coords if hasattr(coord, "mesh") + ] vector_aux_coords = [ - coord for coord in aux_coords - if (id(coord) not in scalar_coord_ids - and coord not in mesh_coords) + coord + for coord in aux_coords + if (id(coord) not in scalar_coord_ids and coord not in mesh_coords) ] vector_derived_coords = [ coord diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py index fe636d1a49..d8ce403f5d 100644 --- a/lib/iris/fileformats/netcdf.py +++ b/lib/iris/fileformats/netcdf.py @@ -2112,13 +2112,18 @@ def _create_mesh(self, cube, mesh): last_dim = f"{cf_mesh_name}_{loc_from}_N_{loc_to}s" # Create if it does not already exist. if last_dim not in self._dataset.dimensions: - length = conn.shape[1] + length = conn.shape[1 - conn.src_dim] self._dataset.createDimension(last_dim, length) # Create variable. # NOTE: for connectivities *with missing points*, this will use a - # fixed standard fill-value of -1, and add a _FillValue property. - conn_dims = (mesh_dims[loc_from], last_dim) + # fixed standard fill-value of -1. In that case, we also add a + # mesh property '_FillValue', below. + loc_dim_name = mesh_dims[loc_from] + conn_dims = (loc_dim_name, last_dim) + if conn.src_dim == 1: + # 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 ) @@ -2128,12 +2133,19 @@ def _create_mesh(self, cube, mesh): _setncattr(cf_conn_var, "start_index", conn.start_index) # If content was masked, also add a "_FillValue" property. - # N.B. for now, we are using -1 as a universal 'safe' value. + # N.B. for now at least, use -1 as a universal 'safe' value. if np.ma.is_masked(conn.indices): - _setncattr(cf_var, "_FillValue", -1) + _setncattr(cf_mesh_var, "_FillValue", -1) # Record the connectivity on the parent mesh var. _setncattr(cf_mesh_var, cf_conn_attr_name, cf_conn_name) + # If the connectivity had the 'alternate' dimension order, add the + # relevant dimension property + if conn.src_dim == 1: + loc_dim_attr = f"{loc_from}_dimension" + # Should only get here once. + assert loc_dim_attr not in cf_mesh_var.ncattrs() + _setncattr(cf_mesh_var, loc_dim_attr, loc_dim_name) return cf_mesh_name diff --git a/lib/iris/tests/integration/experimental/test_ugrid_save.py b/lib/iris/tests/integration/experimental/test_ugrid_save.py index 13ef3956fd..3726d16f4c 100644 --- a/lib/iris/tests/integration/experimental/test_ugrid_save.py +++ b/lib/iris/tests/integration/experimental/test_ugrid_save.py @@ -87,7 +87,7 @@ def test_example_roundtrips(self): if "ex4" in ex_name: # Discard the extra formula terms component cubes # Saving these does not do what you expect - orig_cubes = orig_cubes.extract('datavar') + orig_cubes = orig_cubes.extract("datavar") # Save-and-load-back to compare the Iris saved result. resave_ncfile_path = str(self.temp_dir / f"{ex_name}_resaved.nc") 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 7dcd8b0f9c..55f140362d 100644 --- a/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py +++ b/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py @@ -22,8 +22,6 @@ from iris.cube import Cube, CubeList from iris.experimental.ugrid import Connectivity, Mesh from iris.tests.stock import realistic_4d -from iris.util import new_axis, demote_dim_coord_to_aux_coord - XY_LOCS = ("x", "y") XY_NAMES = ("longitude", "latitude") @@ -115,9 +113,6 @@ def apply_xyargs(coords, xyargs): ) applyargs(mesh, mesh_kwargs) - # N.B. *nasty!* attach an extra convenience property for handling purposes. - mesh._mesh_dims = mesh_dims - return mesh @@ -147,12 +142,29 @@ def make_mesh(basic=True, **kwargs): return mesh +def mesh_location_size(mesh, location): + """Get the size of a location-dimension from a mesh.""" + if location == "node": + # Use a node coordinate (which always exists). + node_coord = mesh.node_coords[0] + result = node_coord.shape[0] + else: + # Use a _node_connectivity, if any. + conn_name = f"{location}_node_connectivity" + conn = getattr(mesh, conn_name, None) + if conn is None: + result = 0 + else: + result = conn.shape[conn.src_dim] + return result + + # Pre-create a simple "standard" test mesh for multiple uses _DEFAULT_MESH = make_mesh() def make_cube(mesh=_DEFAULT_MESH, location="face", **kwargs): - dim = mesh._mesh_dims[location] + dim = mesh_location_size(mesh, location) cube = Cube(np.zeros(dim, np.float32)) for meshco in mesh.to_MeshCoords(location): cube.add_aux_coord(meshco, (0,)) @@ -162,13 +174,13 @@ def make_cube(mesh=_DEFAULT_MESH, location="face", **kwargs): def add_height_dim(cube): - # Add an extra inital dimension to the cube. + # Add an extra inital 'height' dimension onto a cube. cube = cube.copy() # Avoid trashing the input cube. - cube.add_aux_coord(AuxCoord([0.], standard_name='height', units='m')) + cube.add_aux_coord(AuxCoord([0.0], standard_name="height", units="m")) # Make three copies with different heights cubes = [cube.copy() for _ in range(3)] for i_cube, cube in enumerate(cubes): - cube.coord('height').points = [i_cube] + cube.coord("height").points = [i_cube] # Merge to create an additional 'height' dimension. cube = CubeList(cubes).merge_cube() return cube @@ -205,7 +217,7 @@ def vars_w_props(varsdict, **kwargs): """ Subset a vars dict, {name:props}, returning only those where each =, defined by the given keywords. - Except that '="*"' means that an attribute '' merely _exists_. + Except that '="*"' means that '' merely _exists_, with any value. """ @@ -333,7 +345,7 @@ def test_basic_mesh(self): ) # There are no edges. - self.assertNotIn('edge_node_connectivity', mesh_props) + self.assertNotIn("edge_node_connectivity", mesh_props) self.assertEqual( len(vars_w_props(vars, cf_role="edge_node_connectivity")), 0 ) @@ -341,18 +353,22 @@ def test_basic_mesh(self): # The dims are precisely (nodes, faces, nodes-per-face), in that order. self.assertEqual( list(dims.keys()), - ['Mesh2d_nodes', 'Mesh2d_faces', 'Mesh2d_face_N_nodes']) + ["Mesh2d_nodes", "Mesh2d_faces", "Mesh2d_face_N_nodes"], + ) # The variables are (mesh, 2*node-coords, 2*face-coords, face-nodes, data), # in that order self.assertEqual( list(vars.keys()), - ['Mesh2d', - 'node_x', 'node_y', - 'face_x', 'face_y', - 'mesh2d_faces', - 'unknown', - ] + [ + "Mesh2d", + "node_x", + "node_y", + "face_x", + "face_y", + "mesh2d_faces", + "unknown", + ], ) # For definiteness, also check against a full CDL snapshot @@ -411,17 +427,21 @@ def test_multi_cubes_identical_meshes(self): # there are exactly 2 meshes in the file mesh_names = vars_meshnames(vars) - self.assertEqual(sorted(mesh_names), ['Mesh2d', 'Mesh2d_0']) + self.assertEqual(sorted(mesh_names), ["Mesh2d", "Mesh2d_0"]) # they use different dimensions - self.assertEqual(vars_meshdim(vars, 'node', mesh_name='Mesh2d'), - 'Mesh2d_nodes') - self.assertEqual(vars_meshdim(vars, 'face', mesh_name='Mesh2d'), - 'Mesh2d_faces') - 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') + self.assertEqual( + vars_meshdim(vars, "node", mesh_name="Mesh2d"), "Mesh2d_nodes" + ) + self.assertEqual( + vars_meshdim(vars, "face", mesh_name="Mesh2d"), "Mesh2d_faces" + ) + 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" + ) # there are exactly two data-variables with a 'mesh' property mesh_datavars = vars_w_props(vars, mesh="*") @@ -429,15 +449,14 @@ def test_multi_cubes_identical_meshes(self): # the data variables reference the two separate meshes a_props, b_props = vars["a"], vars["b"] - self.assertEqual(a_props['mesh'], 'Mesh2d') - self.assertEqual(a_props['location'], 'face') - self.assertEqual(b_props['mesh'], 'Mesh2d_0') - self.assertEqual(b_props['location'], 'face') + self.assertEqual(a_props["mesh"], "Mesh2d") + self.assertEqual(a_props["location"], "face") + self.assertEqual(b_props["mesh"], "Mesh2d_0") + self.assertEqual(b_props["location"], "face") # the data variables map the appropriate node dimensions - self.assertEqual(a_props["_DIMS"], ['Mesh2d_faces']) - self.assertEqual(b_props["_DIMS"], ['Mesh2d_faces_0']) - + self.assertEqual(a_props["_DIMS"], ["Mesh2d_faces"]) + self.assertEqual(b_props["_DIMS"], ["Mesh2d_faces_0"]) def test_multi_cubes_different_mesh(self): # Check that we can correctly distinguish 2 different meshes. @@ -458,11 +477,11 @@ def test_multi_cubes_different_mesh(self): # the main variables reference the respective meshes, and 'face' location a_props, b_props = vars["a"], vars["b"] - mesh_a, loc_a = a_props['mesh'], a_props['location'] - mesh_b, loc_b = b_props['mesh'], b_props['location'] + mesh_a, loc_a = a_props["mesh"], a_props["location"] + mesh_b, loc_b = b_props["mesh"], b_props["location"] self.assertNotEqual(mesh_a, mesh_b) - self.assertEqual(loc_a, 'face') - self.assertEqual(loc_b, 'face') + self.assertEqual(loc_a, "face") + self.assertEqual(loc_b, "face") def test_nonmesh_dim(self): # Check where the data variable has a 'normal' dim and a mesh dim. @@ -474,16 +493,16 @@ def test_nonmesh_dim(self): # have just 1 mesh, including a face and node coordinates. (mesh_name,) = vars_meshnames(vars) - face_dim = vars_meshdim(vars, 'face', mesh_name) - _ = vars_meshdim(vars, 'node', mesh_name) + face_dim = vars_meshdim(vars, "face", mesh_name) + _ = vars_meshdim(vars, "node", mesh_name) # have just 1 data-variable - ((data_name, data_props),) = vars_w_props(vars, mesh='*').items() + ((data_name, data_props),) = vars_w_props(vars, mesh="*").items() # data maps to the height + mesh dims - self.assertEqual(data_props['_DIMS'], ['height', face_dim]) - self.assertEqual(data_props['mesh'], mesh_name) - self.assertEqual(data_props['location'], 'face') + self.assertEqual(data_props["_DIMS"], ["height", face_dim]) + self.assertEqual(data_props["mesh"], mesh_name) + self.assertEqual(data_props["location"], "face") def test_nonmesh_hybrid_dim(self): # Check a case with a hybrid non-mesh dimension @@ -491,8 +510,12 @@ def test_nonmesh_hybrid_dim(self): # Strip off the time and longtude dims, to make it simpler. cube = cube[0, ..., 0] # Remove all the unwanted coords (also loses the coord-system) - lose_coords = ('time', 'forecast_period', - 'grid_longitude', 'grid_latitude') + lose_coords = ( + "time", + "forecast_period", + "grid_longitude", + "grid_latitude", + ) for coord in lose_coords: cube.remove_coord(coord) @@ -501,10 +524,10 @@ def test_nonmesh_hybrid_dim(self): n_places = cube.shape[i_horizontal_dim] mesh = make_mesh( n_faces=n_places, - n_nodes=30 # arbitrary + unrealistic, but doesn't actually matter + n_nodes=30, # arbitrary + unrealistic, but doesn't actually matter ) # Attach the mesh by adding MeshCoords - for coord in mesh.to_MeshCoords('face'): + for coord in mesh.to_MeshCoords("face"): cube.add_aux_coord(coord, (i_horizontal_dim,)) # Save and snapshot the result @@ -512,34 +535,40 @@ def test_nonmesh_hybrid_dim(self): # have just 1 mesh, including face and node coordinates. (mesh_name,) = vars_meshnames(vars) - face_dim = vars_meshdim(vars, 'face', mesh_name) - _ = vars_meshdim(vars, 'node', mesh_name) + face_dim = vars_meshdim(vars, "face", mesh_name) + _ = vars_meshdim(vars, "node", mesh_name) # have hybrid vertical dimension, with all the usual term variables. - self.assertIn('model_level_number', dims) - vert_vars = list(vars_w_dims(vars, ['model_level_number']).keys()) + self.assertIn("model_level_number", dims) + vert_vars = list(vars_w_dims(vars, ["model_level_number"]).keys()) # The list of file variables mapping the vertical dimensio: # = the data-var, plus all the height terms self.assertEqual( vert_vars, - ['air_potential_temperature', 'model_level_number', - 'level_height', 'level_height_bnds', 'sigma', 'sigma_bnds'] + [ + "air_potential_temperature", + "model_level_number", + "level_height", + "level_height_bnds", + "sigma", + "sigma_bnds", + ], ) # have just 1 data-variable, which maps to hybrid-height and mesh dims - ((data_name, data_props),) = vars_w_props(vars, mesh='*').items() - self.assertEqual(data_props['_DIMS'], ['model_level_number', face_dim]) - self.assertEqual(data_props['mesh'], mesh_name) - self.assertEqual(data_props['location'], 'face') + ((data_name, data_props),) = vars_w_props(vars, mesh="*").items() + self.assertEqual(data_props["_DIMS"], ["model_level_number", face_dim]) + self.assertEqual(data_props["mesh"], mesh_name) + self.assertEqual(data_props["location"], "face") - def test_alternate_dim_order(self): + def test_alternate_cube_dim_order(self): # A cube transposed from the 'usual' order # Should work much the same as the "basic" case. - cube_1 = make_cube(var_name='a') + cube_1 = make_cube(var_name="a") cube_1 = add_height_dim(cube_1) cube_2 = cube_1.copy() - cube_2.var_name = 'b' + cube_2.var_name = "b" cube_2.transpose() # Save and snapshot the result @@ -549,42 +578,63 @@ def test_alternate_dim_order(self): (mesh_name,) = vars_meshnames(vars) # both variables reference the same mesh - v_a, v_b = vars['a'], vars['b'] - self.assertEqual(v_a['mesh'], mesh_name) - self.assertEqual(v_a['location'], 'face') - self.assertEqual(v_b['mesh'], mesh_name) - self.assertEqual(v_b['location'], 'face') + v_a, v_b = vars["a"], vars["b"] + self.assertEqual(v_a["mesh"], mesh_name) + self.assertEqual(v_a["location"], "face") + self.assertEqual(v_b["mesh"], mesh_name) + self.assertEqual(v_b["location"], "face") # Check the var dimensions - self.assertEqual(v_a['_DIMS'], ['height', 'Mesh2d_faces']) - self.assertEqual(v_b['_DIMS'], ['Mesh2d_faces', 'height']) - - # def test_mesh_nonstandard_dims(self): - # # A mesh with connectivities in the 'wrong' order - # # This should record the relevant dimension - # mesh = make_mesh() - # # Get the face-node connectivity - # face_nodes_conn = mesh.face_node_connectivity - # # Transpose it : N.B. this sets src_dim=1, asit should be. - # facenodes_nodeface_conn = face_nodes_conn.transpose() - # # Make a new mesh, based on this. - # mesh2 = Mesh( - # topology_dimension=mesh.topology_dimension, - # node_coords_and_axes=zip(mesh.node_coords, XY_LOCS), - # face_coords_and_axes=zip(mesh.face_coords, XY_LOCS), - # connectivities=[facenodes_nodeface_conn], - # ) - # # N.B. replicate our nasty extra testing property (so we can make_cube) - # mesh2._mesh_dims = mesh._mesh_dims - # - # # Build a cube on the modified mesh - # cube = make_cube(mesh=mesh2) - # - # # Save and snapshot the result - # dims, vars = self.check_save(cube) - # - # t_dbg = 0 - # pass + self.assertEqual(v_a["_DIMS"], ["height", "Mesh2d_faces"]) + self.assertEqual(v_b["_DIMS"], ["Mesh2d_faces", "height"]) + + def test_alternate_connectivity_dim_order(self): + # A mesh with some connectivities in the 'other' order. + # This should also create a property with the dimension name + mesh = make_mesh(n_edges=7) + # Get the face-node and edge-node connectivities + face_nodes_conn = mesh.face_node_connectivity + edge_nodes_conn = mesh.edge_node_connectivity + # Transpose them : N.B. this sets src_dim=1, as it should be. + nodesfirst_faces_conn = face_nodes_conn.transpose() + nodesfirst_edges_conn = edge_nodes_conn.transpose() + # Make a new mesh with both face and edge connectivities 'transposed'. + mesh2 = Mesh( + topology_dimension=mesh.topology_dimension, + node_coords_and_axes=zip(mesh.node_coords, XY_LOCS), + face_coords_and_axes=zip(mesh.face_coords, XY_LOCS), + connectivities=[nodesfirst_faces_conn, nodesfirst_edges_conn], + ) + + # Build a cube on the modified mesh + cube = make_cube(mesh=mesh2) + + # Save and snapshot the result + dims, vars = self.check_save(cube) + + # Check shape and dimensions of the associated connectivity variables. + (mesh_name,) = vars_meshnames(vars) + mesh_props = vars[mesh_name] + faceconn_name = mesh_props["face_node_connectivity"] + edgeconn_name = mesh_props["edge_node_connectivity"] + faceconn_props = vars[faceconn_name] + edgeconn_props = vars[edgeconn_name] + self.assertEqual( + faceconn_props["_DIMS"], ["Mesh_2d_face_N_nodes", "Mesh2d_face"] + ) + self.assertEqual( + edgeconn_props["_DIMS"], ["Mesh_2d_edge_N_nodes", "Mesh2d_edge"] + ) + + # Check the dimension lengths are also as expected + self.assertEqual(dims["Mesh2d_face"], 2) + self.assertEqual(dims["Mesh_2d_face_N_nodes"], 4) + self.assertEqual(dims["Mesh2d_edge"], 7) + self.assertEqual(dims["Mesh_2d_edge_N_nodes"], 2) + + # the mesh has extra location-dimension properties + self.assertEqual(mesh_props["face_dimension"], "Mesh2d_face") + self.assertEqual(mesh_props["edge_dimension"], "Mesh2d_edge") if __name__ == "__main__": From 3892c12fa7592387abf7c2779d066c49df5ea68e Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 14 Sep 2021 11:40:43 +0100 Subject: [PATCH 23/26] Remove debug. --- lib/iris/experimental/ugrid/__init__.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/lib/iris/experimental/ugrid/__init__.py b/lib/iris/experimental/ugrid/__init__.py index 91e5359cc0..e31de8dd59 100644 --- a/lib/iris/experimental/ugrid/__init__.py +++ b/lib/iris/experimental/ugrid/__init__.py @@ -3977,8 +3977,6 @@ def _build_mesh_coords(mesh, cf_var): } mesh_dim_name = locations_dimensions[cf_var.location] # (Only expecting 1 mesh dimension per cf_var). - if mesh_dim_name not in cf_var.dimensions: - oops = True mesh_dim = cf_var.dimensions.index(mesh_dim_name) mesh_coords = mesh.to_MeshCoords(location=cf_var.location) From 0d94eb9a88ae8ce4721995476d3b3b93407e1e2b Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 14 Sep 2021 11:48:00 +0100 Subject: [PATCH 24/26] Move complex additional actions from _create_mesh into _add_mesh. --- lib/iris/fileformats/netcdf.py | 139 +++++++++++++++++---------------- 1 file changed, 71 insertions(+), 68 deletions(-) diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py index d8ce403f5d..30d8531f08 100644 --- a/lib/iris/fileformats/netcdf.py +++ b/lib/iris/fileformats/netcdf.py @@ -1411,6 +1411,76 @@ def _add_mesh(self, cube): cf_mesh_name = self._create_mesh(cube, mesh) self._name_coord_map.append(cf_mesh_name, mesh) + cf_mesh_var = self._dataset.variables[cf_mesh_name] + + # Get the mesh-element dim names. + mesh_dims = self._mesh_dims[mesh] + + # Add all the element coordinate variables. + for location in MESH_LOCATIONS: + coords_meshobj_attr = f"{location}_coords" + coords_file_attr = f"{location}_coordinates" + mesh_coords = getattr(mesh, coords_meshobj_attr, None) + if mesh_coords: + coord_names = [] + for coord in mesh_coords: + if coord is None: + continue # an awkward thing that mesh.coords does + coord_name = self._create_generic_cf_array_var( + cube, [], coord, element_dims=(mesh_dims[location],) + ) + coord_names.append(coord_name) + # Record the coordinates (if any) on the mesh variable. + if coord_names: + coord_names = " ".join(coord_names) + _setncattr(cf_mesh_var, coords_file_attr, coord_names) + + # Add all the connectivity variables. + # pre-fetch the set + ignore "None"s -- looks like a bug ? + conns = [conn for conn in mesh.all_connectivities if conn is not None] + for conn in conns: + # Get the connectivity role, = "{loc1}_{loc2}_connectivity". + cf_conn_attr_name = conn.cf_role + loc_from, loc_to, _ = cf_conn_attr_name.split("_") + # Construct a trailing dimension name. + last_dim = f"{cf_mesh_name}_{loc_from}_N_{loc_to}s" + # Create if it does not already exist. + if last_dim not in self._dataset.dimensions: + length = conn.shape[1 - conn.src_dim] + self._dataset.createDimension(last_dim, length) + + # Create variable. + # NOTE: for connectivities *with missing points*, this will use a + # fixed standard fill-value of -1. In that case, we also add a + # mesh property '_FillValue', below. + loc_dim_name = mesh_dims[loc_from] + conn_dims = (loc_dim_name, last_dim) + if conn.src_dim == 1: + # 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 + ) + # Add essential attributes to the Connectivity variable. + cf_conn_var = self._dataset.variables[cf_conn_name] + _setncattr(cf_conn_var, "cf_role", cf_conn_attr_name) + _setncattr(cf_conn_var, "start_index", conn.start_index) + + # If content was masked, also add a "_FillValue" property. + # N.B. for now at least, use -1 as a universal 'safe' value. + if np.ma.is_masked(conn.indices): + _setncattr(cf_mesh_var, "_FillValue", -1) + + # Record the connectivity on the parent mesh var. + _setncattr(cf_mesh_var, cf_conn_attr_name, cf_conn_name) + # If the connectivity had the 'alternate' dimension order, add the + # relevant dimension property + if conn.src_dim == 1: + loc_dim_attr = f"{loc_from}_dimension" + # Should only get here once. + assert loc_dim_attr not in cf_mesh_var.ncattrs() + _setncattr(cf_mesh_var, loc_dim_attr, loc_dim_name) + return cf_mesh_name def _add_inner_related_vars( @@ -2044,7 +2114,7 @@ def _get_mesh_variable_name(self, mesh): def _create_mesh(self, cube, mesh): """ - Create all the variables associated with a mesh in the netCDF dataset. + Create a mesh variable in the netCDF dataset. Args: @@ -2080,73 +2150,6 @@ def _create_mesh(self, cube, mesh): # Add 'standard' names + units atributes self._set_cf_var_attributes(cf_mesh_var, mesh) - # Get the mesh-element dim names. - mesh_dims = self._mesh_dims[mesh] - # Add the element coordinate variables. - for location in MESH_LOCATIONS: - coords_meshobj_attr = f"{location}_coords" - coords_file_attr = f"{location}_coordinates" - mesh_coords = getattr(mesh, coords_meshobj_attr, None) - if mesh_coords: - coord_names = [] - for coord in mesh_coords: - if coord is None: - continue # an awkward thing that mesh.coords does - coord_name = self._create_generic_cf_array_var( - cube, [], coord, element_dims=(mesh_dims[location],) - ) - coord_names.append(coord_name) - # Record the coordinates (if any) on the mesh variable. - if coord_names: - coord_names = " ".join(coord_names) - _setncattr(cf_mesh_var, coords_file_attr, coord_names) - - # Add all the connectivity variables. - # pre-fetch the set + ignore "None"s -- looks like a bug ? - conns = [conn for conn in mesh.all_connectivities if conn is not None] - for conn in conns: - # Get the connectivity role, = "{loc1}_{loc2}_connectivity". - cf_conn_attr_name = conn.cf_role - loc_from, loc_to, _ = cf_conn_attr_name.split("_") - # Construct a trailing dimension name. - last_dim = f"{cf_mesh_name}_{loc_from}_N_{loc_to}s" - # Create if it does not already exist. - if last_dim not in self._dataset.dimensions: - length = conn.shape[1 - conn.src_dim] - self._dataset.createDimension(last_dim, length) - - # Create variable. - # NOTE: for connectivities *with missing points*, this will use a - # fixed standard fill-value of -1. In that case, we also add a - # mesh property '_FillValue', below. - loc_dim_name = mesh_dims[loc_from] - conn_dims = (loc_dim_name, last_dim) - if conn.src_dim == 1: - # 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 - ) - # Add essential attributes to the Connectivity variable. - cf_conn_var = self._dataset.variables[cf_conn_name] - _setncattr(cf_conn_var, "cf_role", cf_conn_attr_name) - _setncattr(cf_conn_var, "start_index", conn.start_index) - - # If content was masked, also add a "_FillValue" property. - # N.B. for now at least, use -1 as a universal 'safe' value. - if np.ma.is_masked(conn.indices): - _setncattr(cf_mesh_var, "_FillValue", -1) - - # Record the connectivity on the parent mesh var. - _setncattr(cf_mesh_var, cf_conn_attr_name, cf_conn_name) - # If the connectivity had the 'alternate' dimension order, add the - # relevant dimension property - if conn.src_dim == 1: - loc_dim_attr = f"{loc_from}_dimension" - # Should only get here once. - assert loc_dim_attr not in cf_mesh_var.ncattrs() - _setncattr(cf_mesh_var, loc_dim_attr, loc_dim_name) - return cf_mesh_name def _set_cf_var_attributes(self, cf_var, element): From 9d189404a68da91fdd374e4d96ff914c4cc4f189 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 14 Sep 2021 14:56:49 +0100 Subject: [PATCH 25/26] Support+test connectivities with missing points. --- lib/iris/fileformats/netcdf.py | 18 +++--- .../fileformats/netcdf/test_Saver__ugrid.py | 57 +++++++++++++++++++ 2 files changed, 67 insertions(+), 8 deletions(-) diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py index 30d8531f08..ae8144ab63 100644 --- a/lib/iris/fileformats/netcdf.py +++ b/lib/iris/fileformats/netcdf.py @@ -1466,11 +1466,6 @@ def _add_mesh(self, cube): _setncattr(cf_conn_var, "cf_role", cf_conn_attr_name) _setncattr(cf_conn_var, "start_index", conn.start_index) - # If content was masked, also add a "_FillValue" property. - # N.B. for now at least, use -1 as a universal 'safe' value. - if np.ma.is_masked(conn.indices): - _setncattr(cf_mesh_var, "_FillValue", -1) - # Record the connectivity on the parent mesh var. _setncattr(cf_mesh_var, cf_conn_attr_name, cf_conn_name) # If the connectivity had the 'alternate' dimension order, add the @@ -2282,8 +2277,12 @@ def _create_generic_cf_array_var( data = self._ensure_valid_dtype(data, element_type, element) if fill_value is not None: - # Use a specific fill-value in place of the netcdf default. - data = np.ma.filled(data, fill_value) + if np.ma.is_masked(data): + # Use a specific fill-value in place of the netcdf default. + data = np.ma.filled(data, fill_value) + else: + # Create variable without a (non-standard) fill_value + fill_value = None # Check if this is a dim-coord. is_dimcoord = element in cube.dim_coords @@ -2306,7 +2305,10 @@ def _create_generic_cf_array_var( # Create the CF-netCDF variable. cf_var = self._dataset.createVariable( - cf_name, data.dtype.newbyteorder("="), element_dims + cf_name, + data.dtype.newbyteorder("="), + element_dims, + fill_value=fill_value, ) # Add the axis attribute for spatio-temporal CF-netCDF coordinates. 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 55f140362d..9cd00a2676 100644 --- a/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py +++ b/lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py @@ -636,6 +636,63 @@ def test_alternate_connectivity_dim_order(self): self.assertEqual(mesh_props["face_dimension"], "Mesh2d_face") self.assertEqual(mesh_props["edge_dimension"], "Mesh2d_edge") + def test_nonuniform_connectivity(self): + # Check handling of connecitivities with missing points. + n_faces = 7 + mesh = make_mesh(n_faces=n_faces) + + # In this case, add on a partial face-face connectivity. + # construct a vaguely plausible face-face index array + indices = np.ma.arange(n_faces * 4).reshape((7, 4)) + indices = indices % 7 + # include some missing points -- i.e. not all faces have 4 neighbours + indices[(2, (2, 3))] = np.ma.masked + indices[(3, (0, 2))] = np.ma.masked + indices[6, :] = np.ma.masked + + conn = Connectivity( + indices, + cf_role="face_face_connectivity", + ) + mesh.add_connectivities(conn) + cube = make_cube(mesh=mesh) + + # Save and snapshot the result + dims, vars = self.check_save(cube) + + # Check that the mesh saved with the additional connectivity + (mesh_name,) = vars_meshnames(vars) + mesh_props = vars[mesh_name] + self.assertIn("face_face_connectivity", mesh_props) + ff_conn_name = mesh_props["face_face_connectivity"] + + # check that the connectivity has the corrects dims and fill-property + ff_props = vars[ff_conn_name] + self.assertEqual( + ff_props["_DIMS"], ["Mesh2d_faces", "Mesh2d_face_N_faces"] + ) + self.assertIn("_FillValue", ff_props) + self.assertEqual(ff_props["_FillValue"], -1) + + # Check that a 'normal' connectivity does *not* have a _FillValue + fn_conn_name = mesh_props["face_node_connectivity"] + fn_props = vars[fn_conn_name] + self.assertNotIn("_FillValue", fn_props) + + # For what it's worth, *also* check the actual data array in the file + ds = nc.Dataset(self.tempfile_path) + try: + conn_var = ds.variables[ff_conn_name] + data = conn_var[:] + finally: + ds.close() + self.assertIsInstance(data, np.ma.MaskedArray) + self.assertEqual(data.fill_value, -1) + # Compare raw values stored, to indices with -1 at missing points + raw_data = data.data + filled_indices = indices.filled(-1) + self.assertArrayEqual(raw_data, filled_indices) + if __name__ == "__main__": tests.main() From d978f152d6cd4942123775cc8e299f0b873a7ab0 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 14 Sep 2021 23:32:30 +0100 Subject: [PATCH 26/26] Code moved from _create_mesh to _add_mesh needs to be inside conditional. --- lib/iris/fileformats/netcdf.py | 133 +++++++++++++++++---------------- 1 file changed, 69 insertions(+), 64 deletions(-) diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py index ae8144ab63..74215a9ae3 100644 --- a/lib/iris/fileformats/netcdf.py +++ b/lib/iris/fileformats/netcdf.py @@ -1411,70 +1411,75 @@ def _add_mesh(self, cube): cf_mesh_name = self._create_mesh(cube, mesh) self._name_coord_map.append(cf_mesh_name, mesh) - cf_mesh_var = self._dataset.variables[cf_mesh_name] - - # Get the mesh-element dim names. - mesh_dims = self._mesh_dims[mesh] - - # Add all the element coordinate variables. - for location in MESH_LOCATIONS: - coords_meshobj_attr = f"{location}_coords" - coords_file_attr = f"{location}_coordinates" - mesh_coords = getattr(mesh, coords_meshobj_attr, None) - if mesh_coords: - coord_names = [] - for coord in mesh_coords: - if coord is None: - continue # an awkward thing that mesh.coords does - coord_name = self._create_generic_cf_array_var( - cube, [], coord, element_dims=(mesh_dims[location],) - ) - coord_names.append(coord_name) - # Record the coordinates (if any) on the mesh variable. - if coord_names: - coord_names = " ".join(coord_names) - _setncattr(cf_mesh_var, coords_file_attr, coord_names) - - # Add all the connectivity variables. - # pre-fetch the set + ignore "None"s -- looks like a bug ? - conns = [conn for conn in mesh.all_connectivities if conn is not None] - for conn in conns: - # Get the connectivity role, = "{loc1}_{loc2}_connectivity". - cf_conn_attr_name = conn.cf_role - loc_from, loc_to, _ = cf_conn_attr_name.split("_") - # Construct a trailing dimension name. - last_dim = f"{cf_mesh_name}_{loc_from}_N_{loc_to}s" - # Create if it does not already exist. - if last_dim not in self._dataset.dimensions: - length = conn.shape[1 - conn.src_dim] - self._dataset.createDimension(last_dim, length) - - # Create variable. - # NOTE: for connectivities *with missing points*, this will use a - # fixed standard fill-value of -1. In that case, we also add a - # mesh property '_FillValue', below. - loc_dim_name = mesh_dims[loc_from] - conn_dims = (loc_dim_name, last_dim) - if conn.src_dim == 1: - # 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 - ) - # Add essential attributes to the Connectivity variable. - cf_conn_var = self._dataset.variables[cf_conn_name] - _setncattr(cf_conn_var, "cf_role", cf_conn_attr_name) - _setncattr(cf_conn_var, "start_index", conn.start_index) - - # Record the connectivity on the parent mesh var. - _setncattr(cf_mesh_var, cf_conn_attr_name, cf_conn_name) - # If the connectivity had the 'alternate' dimension order, add the - # relevant dimension property - if conn.src_dim == 1: - loc_dim_attr = f"{loc_from}_dimension" - # Should only get here once. - assert loc_dim_attr not in cf_mesh_var.ncattrs() - _setncattr(cf_mesh_var, loc_dim_attr, loc_dim_name) + cf_mesh_var = self._dataset.variables[cf_mesh_name] + + # Get the mesh-element dim names. + mesh_dims = self._mesh_dims[mesh] + + # Add all the element coordinate variables. + for location in MESH_LOCATIONS: + coords_meshobj_attr = f"{location}_coords" + coords_file_attr = f"{location}_coordinates" + mesh_coords = getattr(mesh, coords_meshobj_attr, None) + if mesh_coords: + coord_names = [] + for coord in mesh_coords: + if coord is None: + continue # an awkward thing that mesh.coords does + coord_name = self._create_generic_cf_array_var( + cube, + [], + coord, + element_dims=(mesh_dims[location],), + ) + coord_names.append(coord_name) + # Record the coordinates (if any) on the mesh variable. + if coord_names: + coord_names = " ".join(coord_names) + _setncattr(cf_mesh_var, coords_file_attr, coord_names) + + # Add all the connectivity variables. + # pre-fetch the set + ignore "None"s -- looks like a bug ? + conns = [ + conn for conn in mesh.all_connectivities if conn is not None + ] + for conn in conns: + # Get the connectivity role, = "{loc1}_{loc2}_connectivity". + cf_conn_attr_name = conn.cf_role + loc_from, loc_to, _ = cf_conn_attr_name.split("_") + # Construct a trailing dimension name. + last_dim = f"{cf_mesh_name}_{loc_from}_N_{loc_to}s" + # Create if it does not already exist. + if last_dim not in self._dataset.dimensions: + length = conn.shape[1 - conn.src_dim] + self._dataset.createDimension(last_dim, length) + + # Create variable. + # NOTE: for connectivities *with missing points*, this will use a + # fixed standard fill-value of -1. In that case, we also add a + # mesh property '_FillValue', below. + loc_dim_name = mesh_dims[loc_from] + conn_dims = (loc_dim_name, last_dim) + if conn.src_dim == 1: + # 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 + ) + # Add essential attributes to the Connectivity variable. + cf_conn_var = self._dataset.variables[cf_conn_name] + _setncattr(cf_conn_var, "cf_role", cf_conn_attr_name) + _setncattr(cf_conn_var, "start_index", conn.start_index) + + # Record the connectivity on the parent mesh var. + _setncattr(cf_mesh_var, cf_conn_attr_name, cf_conn_name) + # If the connectivity had the 'alternate' dimension order, add the + # relevant dimension property + if conn.src_dim == 1: + loc_dim_attr = f"{loc_from}_dimension" + # Should only get here once. + assert loc_dim_attr not in cf_mesh_var.ncattrs() + _setncattr(cf_mesh_var, loc_dim_attr, loc_dim_name) return cf_mesh_name