diff --git a/Changelog.rst b/Changelog.rst index 1bcd39f33..159debd63 100644 --- a/Changelog.rst +++ b/Changelog.rst @@ -1,4 +1,13 @@ Version NEXTVERSION +---------------- + +**2026-03-??** + +* Write UGRID datasets with `cfdm.write` + (https://github.com/NCAS-CMS/cfdm/issues/271) + +Version 1.12.4.0 +---------------- **2026-01-??** @@ -9,10 +18,10 @@ Version NEXTVERSION `cfdm.read` (https://github.com/NCAS-CMS/cfdm/issues/355) * New function `cfdm.dataset_flatten` that replaces the deprecated `cfdm.netcdf_flatten` (https://github.com/NCAS-CMS/cfdm/issues/355) -* New optional dependency: ``zarr>=3.1.3`` -* Removed dependency (now optional): ``zarr>=3.0.8`` * Reduce the time taken to import `cfdm` (https://github.com/NCAS-CMS/cfdm/issues/361) +* New optional dependency: ``zarr>=3.1.3`` +* Removed dependency (now optional): ``zarr>=3.0.8`` ---- diff --git a/cfdm/cellconnectivity.py b/cfdm/cellconnectivity.py index bda511bd7..e6a2a4319 100644 --- a/cfdm/cellconnectivity.py +++ b/cfdm/cellconnectivity.py @@ -4,6 +4,7 @@ class CellConnectivity( mixin.NetCDFVariable, + mixin.NetCDFConnectivityDimension, mixin.Topology, mixin.PropertiesData, mixin.Files, diff --git a/cfdm/cfdmimplementation.py b/cfdm/cfdmimplementation.py index 3e3a407d2..c33a1ded1 100644 --- a/cfdm/cfdmimplementation.py +++ b/cfdm/cfdmimplementation.py @@ -419,6 +419,25 @@ def get_bounds_ncvar(self, parent, default=None): return self.nc_get_variable(bounds, default=default) + def get_cell_connectivities(self, parent): + """Return the cell connectivities from a parent. + + .. versionadded:: (cfdm) NEXTVERSION + + :Parameters: + + parent: `Field` or `Domain` + The parent object. + + :Returns: + + `dict` + A dictionary whose values are cell connectivity + objects, keyed by unique identifiers. + + """ + return parent.cell_connectivities(todict=True) + def get_cell_measures(self, field): """Return all of the cell measure constructs of a field. @@ -1012,6 +1031,25 @@ def get_domain_axis_size(self, field, axis): """ return field.domain_axes(todict=True)[axis].get_size() + def get_domain_topologies(self, parent): + """Return the domain topologies from a parent. + + .. versionadded:: (cfdm) NEXTVERSION + + :Parameters: + + parent: `Field` or `Domain` + The parent object. + + :Returns: + + `dict` + A dictionary whose values are domain topology objects, + keyed by unique identifiers. + + """ + return parent.domain_topologies(todict=True) + def get_sample_dimension_position(self, construct): """Returns the position of the compressed data sample dimension. @@ -1983,16 +2021,23 @@ def get_tie_points(self, construct, default=None): return data.source(default) - def initialise_AuxiliaryCoordinate(self): + def initialise_AuxiliaryCoordinate(self, **kwargs): """Return an auxiliary coordinate construct. + :Parameters: + + kwargs: optional + Parameters with which to intialising the object. + + .. versionadded:: (cfdm) NEXTVERSION + :Returns: `AuxiliaryCoordinate` """ cls = self.get_class("AuxiliaryCoordinate") - return cls() + return cls(**kwargs) def initialise_Bounds(self): """Return a bounds component. @@ -3421,26 +3466,6 @@ def set_dependent_tie_points(self, construct, tie_points, dimensions): ) construct.set_data(data) - def set_mesh_id(self, parent, mesh_id): - """Set a mesh identifier. - - .. versionadded:: (cfdm) 1.11.0.0 - - :Parameters: - - parent: construct - The construct on which to set the mesh id - - mesh_id: - The mesh identifier. - - :Returns: - - `None` - - """ - parent.set_mesh_id(mesh_id) - def nc_set_external(self, construct): """Set the external status of a construct. diff --git a/cfdm/data/data.py b/cfdm/data/data.py index b2ae1520e..1dbb1807c 100644 --- a/cfdm/data/data.py +++ b/cfdm/data/data.py @@ -737,12 +737,6 @@ def __getitem__(self, indices): # ------------------------------------------------------------ new._set_dask(dx, clear=self._ALL ^ self._CFA, in_memory=None) - if 0 in new.shape: - raise IndexError( - f"Index [{original_indices}] selects no elements from " - f"data with shape {original_shape}" - ) - # ------------------------------------------------------------ # Get the axis identifiers for the subspace # ------------------------------------------------------------ diff --git a/cfdm/data/mixin/compressedarraymixin.py b/cfdm/data/mixin/compressedarraymixin.py index 5b97fe515..1f20855b7 100644 --- a/cfdm/data/mixin/compressedarraymixin.py +++ b/cfdm/data/mixin/compressedarraymixin.py @@ -5,49 +5,6 @@ class CompressedArrayMixin: """ - def _lock_file_read(self, array): - """Try to return an array that doesn't support concurrent reads. - - .. versionadded:: (cfdm) 1.11.2.0 - - :Parameters: - - array: array_like - The array to process. - - :Returns" - - `dask.array.Array` or array_like - The new `dask` array, or the orginal array if it - couldn't be ascertained how to form the `dask` array. - - """ - try: - return array.to_dask_array() - except AttributeError: - pass - - try: - chunks = array.chunks - except AttributeError: - chunks = "auto" - - try: - array = array.source() - except (ValueError, AttributeError): - pass - - try: - array.get_filename() - except AttributeError: - pass - else: - import dask.array as da - - array = da.from_array(array, chunks=chunks, lock=True) - - return array - def to_dask_array(self, chunks="auto"): """Convert the data to a `dask` array. @@ -87,18 +44,7 @@ def to_dask_array(self, chunks="auto"): context = partial(config.set, scheduler="synchronous") - # If possible, convert the compressed data to a dask array - # that doesn't support concurrent reads. This prevents - # "compute called by compute" failures problems at compute - # time. - # - # TODO: This won't be necessary if this is refactored so that - # the compressed data is part of the same dask graph as - # the compressed subarrays. conformed_data = self.conformed_data() - conformed_data = { - k: self._lock_file_read(v) for k, v in conformed_data.items() - } subarray_kwargs = {**conformed_data, **self.subarray_parameters()} # Get the (cfdm) subarray class diff --git a/cfdm/data/netcdfindexer.py b/cfdm/data/netcdfindexer.py index 4f1efd370..f7d212cd1 100644 --- a/cfdm/data/netcdfindexer.py +++ b/cfdm/data/netcdfindexer.py @@ -249,13 +249,15 @@ def __getitem__(self, index): # E.g. index : (1, np.newaxis, slice(1, 5)) # => index1: (1, slice(1, 5)) # and index2: (slice(None), np.newaxis, slice(None)) - except ValueError: + except (ValueError, TypeError): # Something went wrong, which is indicative of the # variable not supporting the appropriate slicing method # (e.g. `h5netcdf` might have returned "ValueError: Step - # must be >= 1 (got -2)"). Therefore we'll just get the - # entire array as a numpy array, and then try indexing - # that. + # must be >= 1 (got -2)" or "TypeError: Indexing elements + # must be in increasing order"). + # + # Therefore we'll just get the entire array as a numpy + # array, and then try indexing that. data = self._index(Ellipsis) data = self._index(index, data=data) diff --git a/cfdm/data/subarray/cellconnectivitysubarray.py b/cfdm/data/subarray/cellconnectivitysubarray.py index 8dc83b3ff..03b11f41e 100644 --- a/cfdm/data/subarray/cellconnectivitysubarray.py +++ b/cfdm/data/subarray/cellconnectivitysubarray.py @@ -40,6 +40,7 @@ def __getitem__(self, indices): stop += 1 data = self._select_data(check_mask=True) + if np.ma.isMA(data): empty = np.ma.empty else: diff --git a/cfdm/data/subarray/mixin/pointtopology.py b/cfdm/data/subarray/mixin/pointtopology.py index 4567b3ecc..61313f250 100644 --- a/cfdm/data/subarray/mixin/pointtopology.py +++ b/cfdm/data/subarray/mixin/pointtopology.py @@ -25,7 +25,29 @@ def __getitem__(self, indices): from cfdm.functions import integer_dtype start_index = self.start_index - node_connectivity = self._select_data(check_mask=False) + node_connectivity = self._select_data(check_mask=True) + + # ------------------------------------------------------------ + # E.g. For faces, 'node_connectivity' might be (two + # quadrilaterals and one triangle): + # + # [[3 4 2 1 ] + # [5 6 4 3 ] + # [7 2 4 --]] + # + # E.g. For the nine edges of the above faces, + # 'node_connectivity' could be: + # + # [[2 7] + # [4 7] + # [4 2] + # [1 2] + # [3 1] + # [3 4] + # [3 5] + # [6 5] + # [4 6]] + # ------------------------------------------------------------ masked = np.ma.isMA(node_connectivity) @@ -49,8 +71,15 @@ def __getitem__(self, indices): cols_extend = cols.extend u_extend = u.extend + unique_nodes = np.unique(node_connectivity) + if masked: + # Remove the missing value from unique nodes + unique_nodes = unique_nodes[:-1] + + unique_nodes = unique_nodes.tolist() + # WARNING (TODO): This loop is a potential performance bottleneck. - for node in np.unique(node_connectivity).tolist(): + for node in unique_nodes: # Find the collection of all nodes that are joined to this # node via links in the mesh, including this node itself # (which will be at the start of the list). @@ -62,10 +91,11 @@ def __getitem__(self, indices): cols_extend(range(n_nodes)) u_extend(nodes) + del unique_nodes + u = np.array(u, dtype=integer_dtype(largest_node_id)) u = csr_array((u, cols, pointers)) u = u.toarray() - if any(map(isnan, self.shape)): # Store the shape, now that is it known. self._set_component("shape", u.shape, copy=False) @@ -76,6 +106,20 @@ def __getitem__(self, indices): # Mask all zeros u = np.ma.where(u == 0, np.ma.masked, u) + # ------------------------------------------------------------ + # E.g. For either of the face and edge examples above, 'u' + # would now be: + # + # [[1 2 3 -- --] + # [2 1 4 7 --] + # [3 1 4 5 --] + # [4 2 3 6 7] + # [5 3 6 -- --] + # [6 4 5 -- --] + # [7 2 4 -- --]] + # + # ------------------------------------------------------------ + if not start_index: # Subtract 1 to get back to zero-based node identities u -= 1 diff --git a/cfdm/data/subarray/pointtopologyfromfacessubarray.py b/cfdm/data/subarray/pointtopologyfromfacessubarray.py index 7b96cace6..2a5740942 100644 --- a/cfdm/data/subarray/pointtopologyfromfacessubarray.py +++ b/cfdm/data/subarray/pointtopologyfromfacessubarray.py @@ -13,7 +13,8 @@ class PointTopologyFromFacesSubarray(PointTopology, MeshSubarray): """ - def _connected_nodes(self, node, node_connectivity, masked): + @classmethod + def _connected_nodes(self, node, node_connectivity, masked, edges=False): """Return nodes that are joined to *node* by face edges. The input *node* is included at the start of the returned @@ -33,11 +34,32 @@ def _connected_nodes(self, node, node_connectivity, masked): Whether or not *node_connectivity* has masked elements. + edges: `bool`, optional + By default *edges* is False and a flat list of nodes, + including *node* itself at the start, is returned. If + True then a list of edge definitions (i.e. a list of + ordered 2-tuples of nodes) is returned instead. + + .. versionadded:: (cfdm) NEXTVERSION + :Returns: `list` All nodes that are joined to *node*, including *node* - itself at the start. + itself at the start. If *edges* is True then a list of + edge definitions is returned instead. + + **Examples** + + >>> p._connected_nodes(7, nc) + [7, 1, 2, 9] + >>> p._connected_nodes(7, nc, edges=True) + [(1, 7), (2, 7), (7, 9)] + + >>> p._connected_nodes(2, nc) + [2, 7, 8] + >>> p._connected_nodes(2, nc, edges=True) + [(2, 7), (2, 8)] """ if masked: @@ -46,30 +68,42 @@ def _connected_nodes(self, node, node_connectivity, masked): where = np.where # Find the faces that contain this node: - faces = where(node_connectivity == node)[0] + rows, cols = where(node_connectivity == node) nodes = [] nodes_extend = nodes.extend # For each face, find which two of its nodes are neighbours to # 'node'. - for face_nodes in node_connectivity[faces]: + for row, col in zip(node_connectivity[rows], cols): if masked: - face_nodes = face_nodes.compressed() - - face_nodes = face_nodes.tolist() - face_nodes.append(face_nodes[0]) - nodes_extend( - [ - m - for m, n in zip(face_nodes[:-1], face_nodes[1:]) - if n == node - ] - ) - - nodes = list(set(nodes)) - - # Insert 'node' at the front of the list - nodes.insert(0, node) + row = row.compressed() + + row = row.tolist() + + # Get the neighbours of 'node', which is in position 'col' + # in the face. + if not col: + # 'node' is in position 0, so its neighbours are in + # positions -1 and 1 + nodes_extend((row[-1], row[1])) + elif col == len(row) - 1: + # 'node' is in position -1, so its neighbours are in + # positions -2 and 0 + nodes_extend((row[-2], row[0])) + else: + # 'node' is in any other position (col), so its + # neighbours are in positions col-1 and col+1 + nodes_extend((row[col - 1], row[col + 1])) + + nodes = sorted(set(nodes)) + + if edges: + # Return a list of ordered edge definitions + nodes = [(node, n) if node < n else (n, node) for n in nodes] + else: + # Return a flat list of nodes, inserting 'node' at the + # start. + nodes.insert(0, node) return nodes diff --git a/cfdm/data/subsampledarray.py b/cfdm/data/subsampledarray.py index 97c720516..a238419b1 100644 --- a/cfdm/data/subsampledarray.py +++ b/cfdm/data/subsampledarray.py @@ -1293,22 +1293,6 @@ def to_dask_array(self, chunks="auto"): parameters = conformed_data["parameters"] dependent_tie_points = conformed_data["dependent_tie_points"] - # If possible, convert the compressed data, parameters and - # dependent tie points to dask arrays that don't support - # concurrent reads. This prevents "compute called by compute" - # failures problems at compute time. - # - # TODO: This won't be necessary if this is refactored so that - # arrays are part of the same dask graph as the - # compressed subarrays. - compressed_data = self._lock_file_read(compressed_data) - parameters = { - k: self._lock_file_read(v) for k, v in parameters.items() - } - dependent_tie_points = { - k: self._lock_file_read(v) for k, v in dependent_tie_points.items() - } - # Get the (cfdm) subarray class Subarray = self.get_Subarray() subarray_name = Subarray().__class__.__name__ diff --git a/cfdm/domain.py b/cfdm/domain.py index 63ed3ad8c..6974a1ecc 100644 --- a/cfdm/domain.py +++ b/cfdm/domain.py @@ -18,6 +18,7 @@ class Domain( mixin.NetCDFGroupAttributes, mixin.NetCDFComponents, mixin.NetCDFUnreferenced, + mixin.NetCDFMeshVariable, mixin.Properties, mixin.Files, core.Domain, diff --git a/cfdm/domaintopology.py b/cfdm/domaintopology.py index 8d00d567e..7c6eaf277 100644 --- a/cfdm/domaintopology.py +++ b/cfdm/domaintopology.py @@ -1,9 +1,12 @@ +import numpy as np + from . import core, mixin from .decorators import _inplace_enabled, _inplace_enabled_define_and_cleanup class DomainTopology( mixin.NetCDFVariable, + mixin.NetCDFConnectivityDimension, mixin.Topology, mixin.PropertiesData, mixin.Files, @@ -324,6 +327,8 @@ def normalise( .. versionadded:: (cfdm) 1.11.0.0 + .. seealso:: `sort`, `to_edge` + :Parameters: start_index: `int`, optional @@ -400,8 +405,7 @@ def normalise( cell = self.get_cell(None) if cell is None: raise ValueError( - f"Can't normalise {self.__class__.__name__} with unknown " - "cell type" + f"Can't normalise {self!r} with unknown cell type" ) d = _inplace_enabled_define_and_cleanup(self) @@ -425,7 +429,330 @@ def normalise( data, start_index, remove_empty_columns ) else: - raise ValueError(f"Can't normalise: Unknown cell type {cell!r}") + raise ValueError( + f"Can't normalise {self!r}: Unknown cell type {cell!r}" + ) d.set_data(data, copy=False) return d + + def sort(self): + """Sort the domain topology node ids. + + Only edge and point domain topologies can be sorted. + + Sorting is across both array dimensions. In general, dimension + 1 is sorted first, and then dimension 0 is sort by the values + in the first column. + + For an edge domain topology, the second column is also sorted + within each unique first column value. + + For a point dimension topology, the first column is omitted + from the dimension 1 sort (because it contains the node id + definition for each row). + + See the examples for more details. + + .. note:: The purpose of this method is to facilitate the + comparison of normalised domain topologies, to see + if they belong to the same UGRID mesh. The sorted + domain topology will, in general, be inconsistent + with other metadata, such as the node geo-locations + stored as domain cell coordinates or cell + bounds. For this reason, `sort` is not allowed to + occur in-place. + + .. versionadded:: (cfdm) NEXTVERSION + + .. seealso:: `normalise`, `to_edge` + + :Returns: + + `{{class}}` + The sorted domain topology construct. + + **Examples** + + >>> f= {{package}}.example_field(9) + >>> print(f) + Field: northward_wind (ncvar%v) + ------------------------------- + Data : northward_wind(time(2), ncdim%nMesh2_edge(9)) ms-1 + Cell methods : time(2): point (interval: 3600 s) + Dimension coords: time(2) = [2016-01-02 01:00:00, 2016-01-02 11:00:00] gregorian + Auxiliary coords: longitude(ncdim%nMesh2_edge(9)) = [-41.5, ..., -43.0] degrees_east + : latitude(ncdim%nMesh2_edge(9)) = [34.5, ..., 32.0] degrees_north + Topologies : cell:edge(ncdim%nMesh2_edge(9), 2) = [[1, ..., 5]] + Connectivities : connectivity:node(ncdim%nMesh2_edge(9), 6) = [[0, ..., --]] + >>> dt = f.domain_topology() + >>> print(dt.array) + [[1 6] + [3 6] + [3 1] + [0 1] + [2 0] + [2 3] + [2 4] + [5 4] + [3 5]] + >>> print(dt.sort().array) + [[0 1] + [0 2] + [1 3] + [1 6] + [2 3] + [2 4] + [3 5] + [3 6] + [4 5]] + + >>> f= {{package}}.example_field(10) + >>> print(f) + Field: air_pressure (ncvar%pa) + ------------------------------ + Data : air_pressure(time(2), ncdim%nMesh2_node(7)) hPa + Cell methods : time(2): point (interval: 3600 s) + Dimension coords: time(2) = [2016-01-02 01:00:00, 2016-01-02 11:00:00] gregorian + Auxiliary coords: longitude(ncdim%nMesh2_node(7)) = [-45.0, ..., -40.0] degrees_east + : latitude(ncdim%nMesh2_node(7)) = [35.0, ..., 34.0] degrees_north + Topologies : cell:point(ncdim%nMesh2_node(7), 5) = [[0, ..., --]] + >>> dt = f.domain_topology() + >>> dt = dt[::-1] + >>> print(dt.array) + [[6 3 1 -- --] + [5 4 3 -- --] + [4 2 5 -- --] + [3 2 1 5 6] + [2 0 3 4 --] + [1 6 0 3 --] + [0 1 2 -- --]] + >>> print(dt.sort().array) + [[0 1 2 -- --] + [1 0 3 6 --] + [2 0 3 4 --] + [3 1 2 5 6] + [4 2 5 -- --] + [5 3 4 -- --] + [6 1 3 -- --]] + + """ + cell = self.get_cell(None) + if cell not in ("edge", "point"): + raise ValueError(f"Can't sort {self!r} with {cell} cells") + + data = self.array + + if cell == "edge": + # Sort within each row + data.sort(axis=1) + # Sort over rows + data = sorted(data.tolist()) + + elif cell == "point": + # Sort within each row from column 1 + data[:, 1:] = np.ma.sort(data[:, 1:], axis=1, endwith=True) + # Sort over rows by value in column 0 + data = data[data[:, 0].argsort()] + + data = self._Data(data, dtype=self.dtype, chunks=self.data.chunks) + + d = self.copy() + d.set_data(data, copy=False) + + return d + + def to_edge(self, sort=False, face_nodes=None): + """Create a new domain topology of edges. + + The edges will defined from the original domain topology, + either as sides of faces, or the links between nodes, or + copied from the existing edges. + + .. versionadded:: (cfdm) NEXTVERSION + + .. seealso:: `normalise`, `sort` + + :Parameters: + + sort: `bool` + If True then sort output edges. This is equivalent to, + but faster than, setting *sort* to False and sorting + the returned `{{class}}` with its `sort` method. + + face_nodes: `None` or sequence of `int`, optional + The unique node ids for 'face' cells. If `None` (the + default) then the node ids will be inferred from the + data, at some computational expense. The order of the + nodes is immaterial. Providing *face_nodes* is an + optimisation in the case that these are known + externally. + + .. warning:: No checks are carried out to ensure that + the given *face_nodes* match the actual + node ids stored in the data, and a + mis-match will result in an exception or, + worse, an incorrect result. + + :Returns: + + `{{class}}` + The domain topology construct of unique edges. + + **Examples** + + >>> f= {{package}}.example_field(8) + >>> print(f) + Field: air_temperature (ncvar%ta) + --------------------------------- + Data : air_temperature(time(2), ncdim%nMesh2_face(3)) K + Cell methods : time(2): point (interval: 3600 s) + Dimension coords: time(2) = [2016-01-02 01:00:00, 2016-01-02 11:00:00] gregorian + Auxiliary coords: longitude(ncdim%nMesh2_face(3)) = [-44.0, -44.0, -42.0] degrees_east + : latitude(ncdim%nMesh2_face(3)) = [34.0, 32.0, 34.0] degrees_north + Topologies : cell:face(ncdim%nMesh2_face(3), 4) = [[2, ..., --]] + Connectivities : connectivity:edge(ncdim%nMesh2_face(3), 5) = [[0, ..., --]] + >>> dt = f[0].domain_topology() + >>> print(dt.array) + [[2 3 1 0] + [4 5 3 2] + [6 1 3 --]] + >>> edge = dt.to_edge() + >>> edge + + >>> print(edge.array) + [[0 1] + [2 4] + [2 3] + [0 2] + [4 5] + [3 6] + [1 6] + [1 3] + [3 5]] + >>> print(dt.to_edge(sort=True).array) + [[0 1] + [0 2] + [1 3] + [1 6] + [2 3] + [2 4] + [3 5] + [3 6] + [4 5]] + + >>> f= {{package}}.example_field(10) + >>> print(f) + Field: air_pressure (ncvar%pa) + ------------------------------ + Data : air_pressure(time(2), ncdim%nMesh2_node(7)) hPa + Cell methods : time(2): point (interval: 3600 s) + Dimension coords: time(2) = [2016-01-02 01:00:00, 2016-01-02 11:00:00] gregorian + Auxiliary coords: longitude(ncdim%nMesh2_node(7)) = [-45.0, ..., -40.0] degrees_east + : latitude(ncdim%nMesh2_node(7)) = [35.0, ..., 34.0] degrees_north + Topologies : cell:point(ncdim%nMesh2_node(7), 5) = [[0, ..., --]] + >>> dt = f[0].domain_topology() + >>> print(dt.array) + [[0 1 2 -- --] + [1 6 0 3 --] + [2 0 3 4 --] + [3 2 1 5 6] + [4 2 5 -- --] + [5 4 3 -- --] + [6 3 1 -- --]] + >>> dt.to_edge() + + >>> print(dt.to_edge(sort=True).array) + [[0 1] + [0 2] + [1 3] + [1 6] + [2 3] + [2 4] + [3 5] + [3 6] + [4 5]] + + """ + edges = [] + edges_extend = edges.extend + + cell = self.get_cell(None) + if face_nodes is not None and cell != "face": + raise ValueError( + "Can't set 'face_nodes' for {self!r} with {cell} cells" + ) + + # Deal with simple "edge" case first + if cell == "edge": + if sort: + edges = self.sort() + else: + edges = self.copy() + + return edges + + # Still here? Then deal with the other cell types. + if cell == "point": + points = self.array + masked = np.ma.is_masked(points) + + # Loop round the nodes, finding the node-pairs that define + # the edges. + for row in points: + if masked: + row = row.compressed() + + node = int(row[0]) + row = row[1:].tolist() + edges_extend( + [(node, n) if node < n else (n, node) for n in row] + ) + + del points, row + + elif cell == "face": + from cfdm.data.subarray import PointTopologyFromFacesSubarray + + connected_nodes = PointTopologyFromFacesSubarray._connected_nodes + + faces = self.array + masked = np.ma.is_masked(faces) + + if face_nodes is None: + # Find the unique node ids + face_nodes = np.unique(faces).tolist() + if masked: + face_nodes = face_nodes[:-1] + + # Loop round the face nodes, finding the node-pairs that + # define the edges. + for n in face_nodes: + edges_extend(connected_nodes(n, faces, masked, edges=True)) + + del faces, face_nodes + + else: + raise NotImplementedError( + f"Can't get edges from {self!r} with {cell} cells" + ) + + # Remove duplicates to get the set of unique edges, because + # every edge currently appears twice in the 'edges' list. + # + # E.g. edge (1, 5) will appear once from processing node 1, + # and once from processing node 5. + edges = set(edges) + + if sort: + edges = sorted(edges) + else: + edges = list(edges) + + edges = self._Data(edges, dtype=self.dtype) + + out = self.copy() + out.set_cell("edge") + out.set_data(edges, copy=False) + + return out diff --git a/cfdm/examplefield.py b/cfdm/examplefield.py index ffcf7f26d..fa4948e14 100644 --- a/cfdm/examplefield.py +++ b/cfdm/examplefield.py @@ -257,8 +257,6 @@ def example_field(n, _implementation=_implementation): Data = _implementation.get_class("Data") - mesh_id = "f51e5aa5e2b0439f9fae4f04e51556f7" - if n == 0: f = Field() @@ -4681,7 +4679,6 @@ def example_field(n, _implementation=_implementation): dtype="f8", ) f.set_data(data) - f.set_mesh_id(mesh_id) # # domain_axis: ncdim%time c = DomainAxis() @@ -4737,7 +4734,7 @@ def example_field(n, _implementation=_implementation): data = Data([-44, -44, -42], units="degrees_east", dtype="f8") c.set_data(data) b = Bounds() - b.nc_set_variable("Mesh2_node_x_bounds") + b.nc_set_variable("Mesh2_face_x_bounds") data = Data( [ [-45, -43, -43, -45], @@ -4763,7 +4760,7 @@ def example_field(n, _implementation=_implementation): data = Data([34, 32, 34], units="degrees_north", dtype="f8") c.set_data(data) b = Bounds() - b.nc_set_variable("Mesh2_node_y_bounds") + b.nc_set_variable("Mesh2_face_y_bounds") data = Data( [ [33, 33, 35, 35], @@ -4785,7 +4782,7 @@ def example_field(n, _implementation=_implementation): c.set_properties({"long_name": "Maps every face to its corner nodes"}) c.nc_set_variable("Mesh2_face_nodes") data = Data( - [[2, 3, 1, 0], [4, 5, 3, 2], [1, 3, 6, -99]], + [[2, 3, 1, 0], [4, 5, 3, 2], [6, 1, 3, -99]], dtype="i4", mask_value=-99, ) @@ -4844,7 +4841,6 @@ def example_field(n, _implementation=_implementation): dtype="f8", ) f.set_data(data) - f.set_mesh_id(mesh_id) # # domain_axis: ncdim%time c = DomainAxis() @@ -4904,7 +4900,7 @@ def example_field(n, _implementation=_implementation): ) c.set_data(data) b = Bounds() - b.nc_set_variable("Mesh2_node_x_bounds") + b.nc_set_variable("Mesh2_edge_x_bounds") data = Data( [ [-43, -40], @@ -4939,7 +4935,7 @@ def example_field(n, _implementation=_implementation): ) c.set_data(data) b = Bounds() - b.nc_set_variable("Mesh2_node_y_bounds") + b.nc_set_variable("Mesh2_edge_y_bounds") data = Data( [ [35, 34], @@ -5050,7 +5046,6 @@ def example_field(n, _implementation=_implementation): dtype="f8", ) f.set_data(data) - f.set_mesh_id(mesh_id) # # domain_axis: ncdim%time c = DomainAxis() @@ -5138,10 +5133,10 @@ def example_field(n, _implementation=_implementation): data = Data( [ [0, 1, 2, -99, -99], - [1, 0, 3, 6, -99], - [2, 0, 4, 3, -99], + [1, 6, 0, 3, -99], + [2, 0, 3, 4, -99], [3, 2, 1, 5, 6], - [4, 5, 2, -99, -99], + [4, 2, 5, -99, -99], [5, 4, 3, -99, -99], [6, 3, 1, -99, -99], ], diff --git a/cfdm/field.py b/cfdm/field.py index 7afe07e5a..00151ff94 100644 --- a/cfdm/field.py +++ b/cfdm/field.py @@ -40,6 +40,7 @@ class Field( mixin.NetCDFGroupAttributes, mixin.NetCDFComponents, mixin.NetCDFUnreferenced, + mixin.NetCDFMeshVariable, mixin.FieldDomain, mixin.PropertiesData, mixin.Files, @@ -272,10 +273,6 @@ def __getitem__(self, indices): if indices is Ellipsis: return new - # Remove a mesh id, on the assumption that the subspaced - # domain will be different to the original. - new.del_mesh_id(None) - data = self.get_data(_fill_value=False) indices = parse_indices(data.shape, indices) @@ -1713,10 +1710,6 @@ def creation_commands( header=header, ) - mesh_id = self.get_mesh_id(None) - if mesh_id is not None: - out.append(f"{name}.set_mesh_id({mesh_id!r})") - nc_global_attributes = self.nc_global_attributes() if nc_global_attributes: out.append("#") diff --git a/cfdm/mixin/__init__.py b/cfdm/mixin/__init__.py index 027bc77e2..ab1285d21 100644 --- a/cfdm/mixin/__init__.py +++ b/cfdm/mixin/__init__.py @@ -13,6 +13,7 @@ from .netcdf import ( NetCDFComponents, + NetCDFConnectivityDimension, NetCDFGlobalAttributes, NetCDFGroupAttributes, NetCDFUnlimitedDimension, @@ -22,6 +23,7 @@ NetCDFGroupsMixin, NetCDFChunks, NetCDFInterpolationSubareaDimension, + NetCDFMeshVariable, NetCDFMixin, NetCDFNodeCoordinateVariable, NetCDFSampleDimension, diff --git a/cfdm/mixin/fielddomain.py b/cfdm/mixin/fielddomain.py index ed5363c57..3e5fb01bb 100644 --- a/cfdm/mixin/fielddomain.py +++ b/cfdm/mixin/fielddomain.py @@ -1,6 +1,7 @@ import logging from ..decorators import _manage_log_level_via_verbosity +from ..functions import _DEPRECATION_ERROR_METHOD logger = logging.getLogger(__name__) @@ -12,39 +13,6 @@ class FieldDomain: """ - def __initialise_from_source(self, source, copy=True): - """Initialise mesh_id information from a source. - - This method is called by - `_Container__parent_initialise_from_source`, which in turn is - called by `cfdm.core.Container.__init__`. - - .. versionadded:: (cfdm) 1.12.2.0 - - :Parameters: - - source: - The object from which to extract the initialisation - information. Typically, but not necessarily, a - `{{class}}` object. - - copy: `bool`, optional - If True (the default) then deep copy the - initialisation information. - - :Returns: - - `None` - - """ - try: - mesh_id = source.get_mesh_id(None) - except AttributeError: - pass - else: - if mesh_id is not None: - self.set_mesh_id(mesh_id) - def _apply_masking_constructs(self): """Apply masking to metadata constructs in-place. @@ -706,7 +674,12 @@ def del_mesh_id(self, default=ValueError()): None """ - return self._del_component("mesh_id", default=default) + _DEPRECATION_ERROR_METHOD( + self, + "del_mesh_id", + version="NEXTVERSION", + removed_at="5.0.0", + ) # pragma: no cover def domain_topology( self, @@ -2315,7 +2288,12 @@ def get_mesh_id(self, default=ValueError()): None """ - return self._get_component("mesh_id", default=default) + _DEPRECATION_ERROR_METHOD( + self, + "get_mesh_id", + version="NEXTVERSION", + removed_at="5.0.0", + ) # pragma: no cover def has_construct(self, *identity, **filter_kwargs): """Whether a unique metadata construct exists. @@ -2434,7 +2412,12 @@ def has_mesh_id(self): None """ - return self._has_component("mesh_id") + _DEPRECATION_ERROR_METHOD( + self, + "has_mesh_id", + version="NEXTVERSION", + removed_at="5.0.0", + ) # pragma: no cover def set_mesh_id(self, mesh_id): """Set a UGRID mesh topology identifier. @@ -2474,4 +2457,9 @@ def set_mesh_id(self, mesh_id): None """ - return self._set_component("mesh_id", mesh_id, copy=False) + _DEPRECATION_ERROR_METHOD( + self, + "set_mesh_id", + version="NEXTVERSION", + removed_at="5.0.0", + ) # pragma: no cover diff --git a/cfdm/mixin/netcdf.py b/cfdm/mixin/netcdf.py index 2da9e9ee1..9c712fa12 100644 --- a/cfdm/mixin/netcdf.py +++ b/cfdm/mixin/netcdf.py @@ -5307,3 +5307,697 @@ def nc_set_dataset_shards(self, shards): ) self._set_netcdf("dataset_shards", shards) + + +class NetCDFMeshVariable(NetCDFMixin, NetCDFGroupsMixin): + """Mixin for accessing the netCDF mesh variable name. + + .. versionadded:: (cfdm) NEXTVERSION + + """ + + def nc_del_mesh_variable(self, default=ValueError()): + """Remove the netCDF mesh variable name. + + .. versionadded:: (cfdm) NEXTVERSION + + .. seealso:: `nc_get_mesh_variable`, `nc_has_mesh_variable`, + `nc_set_mesh_variable` + + :Parameters: + + default: optional + Return the value of the *default* parameter if the + netCDF mesh variable name has not been set. If set to + an `Exception` instance then it will be raised + instead. + + :Returns: + + `str` + The removed netCDF mesh variable name. + + **Examples** + + >>> f.nc_set_mesh_variable('mesh1') + >>> f.nc_has_mesh_variable() + True + >>> f.nc_get_mesh_variable() + 'mesh1' + >>> f.nc_del_mesh_variable() + 'mesh1' + >>> f.nc_has_mesh_variable() + False + >>> print(f.nc_get_mesh_variable(None)) + None + >>> print(f.nc_del_mesh_variable(None)) + None + + """ + return self._nc_del("mesh_variable", default=default) + + def nc_get_mesh_variable(self, default=ValueError()): + """Return the netCDF mesh variable name. + + .. versionadded:: (cfdm) NEXTVERSION + + .. seealso:: `nc_del_mesh_variable`, `nc_has_mesh_variable`, + `nc_set_mesh_variable` + + :Parameters: + + default: optional + Return the value of the *default* parameter if the + netCDF mesh variable name has not been + set. If set to an `Exception` instance then it will be + raised instead. + + :Returns: + + `str` + The netCDF mesh variable name. If unset + then *default* is returned, if provided. + + **Examples** + + >>> f.nc_set_mesh_variable('mesh1') + >>> f.nc_has_mesh_variable() + True + >>> f.nc_get_mesh_variable() + 'mesh1' + >>> f.nc_del_mesh_variable() + 'mesh1' + >>> f.nc_has_mesh_variable() + False + >>> print(f.nc_get_mesh_variable(None)) + None + >>> print(f.nc_del_mesh_variable(None)) + None + + """ + return self._nc_get("mesh_variable", default=default) + + def nc_has_mesh_variable(self): + """Whether the netCDF mesh variable name is set. + + .. versionadded:: (cfdm) NEXTVERSION + + .. seealso:: `nc_get_mesh_variable`, `nc_del_mesh_variable`, + `nc_set_mesh_variable` + + :Returns: + + `bool` + `True` if the netCDF mesh variable name has + been set, otherwise `False`. + + **Examples** + + >>> f.nc_set_mesh_variable('mesh1') + >>> f.nc_has_mesh_variable() + True + >>> f.nc_get_mesh_variable() + 'mesh1' + >>> f.nc_del_mesh_variable() + 'mesh1' + >>> f.nc_has_mesh_variable() + False + >>> print(f.nc_get_mesh_variable(None)) + None + >>> print(f.nc_del_mesh_variable(None)) + None + + """ + return "mesh_variable" in self._get_netcdf() + + def nc_set_mesh_variable(self, value): + """Set the netCDF mesh variable name. + + If there are any ``/`` (slash) characters in the netCDF name + then these act as delimiters for a group hierarchy. By + default, or if the name starts with a ``/`` character and + contains no others, the name is assumed to be in the root + group. + + .. versionadded:: (cfdm) NEXTVERSION + + .. seealso:: `nc_get_mesh_variable`, `nc_has_mesh_variable`, + `nc_del_mesh_variable` + + :Parameters: + + value: `str` + The value for the netCDF mesh variable + name. + + :Returns: + + `None` + + **Examples** + + >>> f.nc_set_mesh_variable('mesh1') + >>> f.nc_has_mesh_variable() + True + >>> f.nc_get_mesh_variable() + 'mesh1' + >>> f.nc_del_mesh_variable() + 'mesh1' + >>> f.nc_has_mesh_variable() + False + >>> print(f.nc_get_mesh_variable(None)) + None + >>> print(f.nc_del_mesh_variable(None)) + None + + """ + return self._nc_set("mesh_variable", value) + + def nc_mesh_variable_groups(self): + """Return the netCDF mesh variable group hierarchy. + + The group hierarchy is defined by the netCDF name. Groups are + delimited by ``/`` (slash) characters in the netCDF name. The + groups are returned, in hierarchical order, as a sequence of + strings. If the name is not set, or contains no ``/`` + characters then an empty sequence is returned, signifying the + root group. + + .. versionadded:: (cfdm) NEXTVERSION + + .. seealso:: `nc_clear_mesh_variable_groups`, + `nc_set_mesh_variable_groups` + + :Returns: + + `tuple` of `str` + The group structure. + + **Examples** + + >>> f.nc_set_mesh_variable('time') + >>> f.nc_mesh_variable_groups() + () + >>> f.nc_set_mesh_variable_groups(['forecast', 'model']) + >>> f.nc_mesh_variable_groups() + ('forecast', 'model') + >>> f.nc_get_mesh_variable() + '/forecast/model/time' + >>> f.nc_clear_mesh_variable_groups() + ('forecast', 'model') + >>> f.nc_get_variable() + 'time' + + >>> f.nc_set_mesh_variable('/forecast/model/time') + >>> f.nc_variable_groups() + ('forecast', 'model') + >>> f.nc_del_mesh_variable('/forecast/model/time') + '/forecast/model/time' + >>> f.nc_mesh_variable_groups() + () + + """ + return self._nc_groups(nc_get=self.nc_get_mesh_variable) + + def nc_set_mesh_variable_groups(self, groups): + """Set the netCDF mesh_variable group hierarchy. + + The group hierarchy is defined by the netCDF name. Groups are + delimited by ``/`` (slash) characters in the netCDF name. The + groups are returned, in hierarchical order, as a sequence of + strings. If the name is not set, or contains no ``/`` + characters then an empty sequence is returned, signifying the + root group. + + An alternative technique for setting the group structure is to + set the netCDF mesh variable name, with + `nc_set_mesh_variable`, with the group structure + delimited by ``/`` characters. + + .. versionadded:: (cfdm) NEXTVERSION + + .. seealso:: `nc_clear_mesh_variable_groups`, + `nc_mesh_variable_groups` + + :Parameters: + + groups: sequence of `str` + The new group structure. + + :Returns: + + `tuple` of `str` + The group structure prior to being reset. + + **Examples** + + >>> f.nc_set_mesh_variable('time') + >>> f.nc_mesh_variable_groups() + () + >>> f.nc_set_mesh_variable_groups(['forecast', 'model']) + >>> f.nc_mesh_variable_groups() + ('forecast', 'model') + >>> f.nc_get_mesh_variable() + '/forecast/model/time' + >>> f.nc_clear_mesh_variable_groups() + ('forecast', 'model') + >>> f.nc_get_variable() + 'time' + + >>> f.nc_set_mesh_variable('/forecast/model/time') + >>> f.nc_variable_groups() + ('forecast', 'model') + >>> f.nc_del_mesh_variable('/forecast/model/time') + '/forecast/model/time' + >>> f.nc_mesh_variable_groups() + () + + """ + return self._nc_set_groups( + groups, + nc_get=self.nc_get_mesh_variable, + nc_set=self.nc_set_mesh_variable, + nc_groups=self.nc_mesh_variable_groups, + ) + + def nc_clear_mesh_variable_groups(self): + """Remove the netCDF mesh variable group hierarchy. + + The group hierarchy is defined by the netCDF name. Groups are + delimited by ``/`` (slash) characters in the netCDF name. The + groups are returned, in hierarchical order, as a sequence of + strings. If the name is not set, or contains no ``/`` + characters then an empty sequence is returned, signifying the + root group. + + An alternative technique for removing the group structure is + to set the netCDF mesh variable name, with + `nc_set_mesh_variable`, with no ``/`` characters. + + .. versionadded:: (cfdm) NEXTVERSION + + .. seealso:: `nc_mesh_variable_groups`, + `nc_set_mesh_variable_groups` + + :Returns: + + `tuple` of `str` + The removed group structure. + + **Examples** + + >>> f.nc_set_mesh_variable('time') + >>> f.nc_mesh_variable_groups() + () + >>> f.nc_set_mesh_variable_groups(['forecast', 'model']) + >>> f.nc_mesh_variable_groups() + ('forecast', 'model') + >>> f.nc_get_mesh_variable() + '/forecast/model/time' + >>> f.nc_clear_mesh_variable_groups() + ('forecast', 'model') + >>> f.nc_get_variable() + 'time' + + >>> f.nc_set_mesh_variable('/forecast/model/time') + >>> f.nc_variable_groups() + ('forecast', 'model') + >>> f.nc_del_mesh_variable('/forecast/model/time') + '/forecast/model/time' + >>> f.nc_mesh_variable_groups() + () + + """ + return self._nc_clear_groups( + nc_get=self.nc_get_mesh_variable, + nc_set=self.nc_set_mesh_variable, + nc_groups=self.nc_mesh_variable_groups, + ) + + +class NetCDFConnectivityDimension(NetCDFMixin, NetCDFGroupsMixin): + """Mixin class for accessing the netCDF connectivity dimension name. + + .. versionadded:: (cfdm) NEXTVERSION + + """ + + def nc_del_connectivity_dimension(self, default=ValueError()): + """Remove the netCDF connectivity dimension name. + + .. versionadded:: (cfdm) NEXTVERSION + + .. seealso:: `nc_get_connectivity_dimension`, + `nc_has_connectivity_dimension`, + `nc_set_connectivity_dimension` + + :Parameters: + + default: optional + Return the value of the *default* parameter if the + netCDF connectivity dimension name has not been + set. If set to an `Exception` instance then it will be + raised instead. + + :Returns: + + `str` + The removed netCDF connectivity dimension name. + + **Examples** + + >>> f.nc_set_connectivity_dimension('time') + >>> f.nc_has_connectivity_dimension() + True + >>> f.nc_get_connectivity_dimension() + 'time' + >>> f.nc_del_connectivity_dimension() + 'time' + >>> f.nc_has_connectivity_dimension() + False + >>> print(f.nc_get_connectivity_dimension(None)) + None + >>> print(f.nc_del_connectivity_dimension(None)) + None + + """ + try: + return self._get_netcdf().pop("connectivity_dimension") + except KeyError: + if default is None: + return default + + return self._default( + default, + f"{self.__class__.__name__} has no netCDF connectivity " + "dimension name", + ) + + def nc_get_connectivity_dimension(self, default=ValueError()): + """Return the netCDF connectivity dimension name. + + .. versionadded:: (cfdm) NEXTVERSION + + .. seealso:: `nc_del_connectivity_dimension`, + `nc_has_connectivity_dimension`, + `nc_set_connectivity_dimension` + + :Parameters: + + default: optional + Return the value of the *default* parameter if the + netCDF connectivity dimension name has not been + set. If set to an `Exception` instance then it will be + raised instead. + + :Returns: + + `str` + The netCDF connectivity dimension name. + + **Examples** + + >>> f.nc_set_connectivity_dimension('time') + >>> f.nc_has_connectivity_dimension() + True + >>> f.nc_get_connectivity_dimension() + 'time' + >>> f.nc_del_connectivity_dimension() + 'time' + >>> f.nc_has_connectivity_dimension() + False + >>> print(f.nc_get_connectivity_dimension(None)) + None + >>> print(f.nc_del_connectivity_dimension(None)) + None + + """ + try: + return self._get_netcdf()["connectivity_dimension"] + except KeyError: + if default is None: + return default + + return self._default( + default, + f"{self.__class__.__name__} has no netCDF connectivity " + "dimension name", + ) + + def nc_has_connectivity_dimension(self): + """Whether the netCDF connectivity dimension name has been set. + + .. versionadded:: (cfdm) NEXTVERSION + + .. seealso:: `nc_del_connectivity_dimension`, + `nc_get_connectivity_dimension`, + `nc_set_connectivity_dimension` + + :Returns: + + `bool` + `True` if the netCDF connectivity dimension name has + been set, otherwise `False`. + + **Examples** + + >>> f.nc_set_connectivity_dimension('time') + >>> f.nc_has_connectivity_dimension() + True + >>> f.nc_get_connectivity_dimension() + 'time' + >>> f.nc_del_connectivity_dimension() + 'time' + >>> f.nc_has_connectivity_dimension() + False + >>> print(f.nc_get_connectivity_dimension(None)) + None + >>> print(f.nc_del_connectivity_dimension(None)) + None + + """ + return "connectivity_dimension" in self._get_netcdf() + + def nc_set_connectivity_dimension(self, value): + """Set the netCDF connectivity dimension name. + + If there are any ``/`` (slash) characters in the netCDF name + then these act as delimiters for a group hierarchy. By + default, or if the name starts with a ``/`` character and + contains no others, the name is assumed to be in the root + group. + + .. versionadded:: (cfdm) NEXTVERSION + + .. seealso:: `nc_del_connectivity_dimension`, + `nc_get_connectivity_dimension`, + `nc_has_connectivity_dimension` + + :Parameters: + + value: `str` + The value for the netCDF connectivity dimension name. + + :Returns: + + `None` + + **Examples** + + >>> f.nc_set_connectivity_dimension('time') + >>> f.nc_has_connectivity_dimension() + True + >>> f.nc_get_connectivity_dimension() + 'time' + >>> f.nc_del_connectivity_dimension() + 'time' + >>> f.nc_has_connectivity_dimension() + False + >>> print(f.nc_get_connectivity_dimension(None)) + None + >>> print(f.nc_del_connectivity_dimension(None)) + None + + """ + if not value or value == "/": + raise ValueError( + f"Invalid netCDF connectivity dimension name: {value!r}" + ) + + if "/" in value: + if not value.startswith("/"): + raise ValueError( + "A netCDF connectivity dimension name with a group " + f"structure must start with a '/'. Got {value!r}" + ) + + if value.count("/") == 1: + value = value[1:] + elif value.endswith("/"): + raise ValueError( + "A netCDF connectivity dimension name with a group " + f"structure can't end with a '/'. Got {value!r}" + ) + + self._set_netcdf("connectivity_dimension", value) + + def nc_connectivity_dimension_groups(self): + """Return the netCDF connectivity dimension group hierarchy. + + The group hierarchy is defined by the netCDF name. Groups are + delimited by ``/`` (slash) characters in the netCDF name. The + groups are returned, in hierarchical order, as a sequence of + strings. If the name is not set, or contains no ``/`` + characters then an empty sequence is returned, signifying the + root group. + + .. versionadded:: (cfdm) 1.8.6 + + .. seealso:: `nc_clear_connectivity_dimension_groups`, + `nc_set_connectivity_dimension_groups` + + :Returns: + + `tuple` of `str` + The group structure. + + **Examples** + + >>> f.nc_set_connectivity_dimension('element') + >>> f.nc_connectivity_dimension_groups() + () + >>> f.nc_set_connectivity_dimension_groups(['forecast', 'model']) + >>> f.nc_connectivity_dimension_groups() + ('forecast', 'model') + >>> f.nc_get_connectivity_dimension() + '/forecast/model/element' + >>> f.nc_clear_connectivity_dimension_groups() + ('forecast', 'model') + >>> f.nc_get_connectivity_dimension() + 'element' + + >>> f.nc_set_connectivity_dimension('/forecast/model/element') + >>> f.nc_connectivity_dimension_groups() + ('forecast', 'model') + >>> f.nc_del_connectivity_dimension('/forecast/model/element') + '/forecast/model/element' + >>> f.nc_connectivity_dimension_groups() + () + + """ + return self._nc_groups(nc_get=self.nc_get_connectivity_dimension) + + def nc_set_connectivity_dimension_groups(self, groups): + """Set the netCDF connectivity dimension group hierarchy. + + The group hierarchy is defined by the netCDF name. Groups are + delimited by ``/`` (slash) characters in the netCDF name. The + groups are returned, in hierarchical order, as a sequence of + strings. If the name is not set, or contains no ``/`` + characters then an empty sequence is returned, signifying the + root group. + + An alternative technique for setting the group structure is to + set the netCDF dimension name, with + `nc_set_connectivity_dimension`, with the group structure + delimited by ``/`` characters. + + .. versionadded:: (cfdm) 1.8.6 + + .. seealso:: `nc_clear_connectivity_dimension_groups`, + `nc_connectivity_dimension_groups` + + :Parameters: + + groups: sequence of `str` + The new group structure. + + :Returns: + + `tuple` of `str` + The group structure prior to being reset. + + **Examples** + + >>> f.nc_set_connectivity_dimension('element') + >>> f.nc_connectivity_dimension_groups() + () + >>> f.nc_set_connectivity_dimension_groups(['forecast', 'model']) + >>> f.nc_connectivity_dimension_groups() + ('forecast', 'model') + >>> f.nc_get_connectivity_dimension() + '/forecast/model/element' + >>> f.nc_clear_connectivity_dimension_groups() + ('forecast', 'model') + >>> f.nc_get_connectivity_dimension() + 'element' + + >>> f.nc_set_connectivity_dimension('/forecast/model/element') + >>> f.nc_connectivity_dimension_groups() + ('forecast', 'model') + >>> f.nc_del_connectivity_dimension('/forecast/model/element') + '/forecast/model/element' + >>> f.nc_connectivity_dimension_groups() + () + + """ + return self._nc_set_groups( + groups, + nc_get=self.nc_get_connectivity_dimension, + nc_set=self.nc_set_connectivity_dimension, + nc_groups=self.nc_connectivity_dimension_groups, + ) + + def nc_clear_connectivity_dimension_groups(self): + """Remove the netCDF connectivity dimension group hierarchy. + + The group hierarchy is defined by the netCDF name. Groups are + delimited by ``/`` (slash) characters in the netCDF name. The + groups are returned, in hierarchical order, as a sequence of + strings. If the name is not set, or contains no ``/`` + characters then an empty sequence is returned, signifying the + root group. + + An alternative technique for removing the group structure is + to set the netCDF dimension name, with + `nc_set_connectivity_dimension`, with no ``/`` characters. + + .. versionadded:: (cfdm) 1.8.6 + + .. seealso:: `nc_connectivity_dimension_groups`, + `nc_set_connectivity_dimension_groups` + + :Returns: + + `tuple` of `str` + The removed group structure. + + **Examples** + + >>> f.nc_set_connectivity_dimension('element') + >>> f.nc_connectivity_dimension_groups() + () + >>> f.nc_set_connectivity_dimension_groups(['forecast', 'model']) + >>> f.nc_connectivity_dimension_groups() + ('forecast', 'model') + >>> f.nc_get_connectivity_dimension() + '/forecast/model/element' + >>> f.nc_clear_connectivity_dimension_groups() + ('forecast', 'model') + >>> f.nc_get_connectivity_dimension() + 'element' + + >>> f.nc_set_connectivity_dimension('/forecast/model/element') + >>> f.nc_connectivity_dimension_groups() + ('forecast', 'model') + >>> f.nc_del_connectivity_dimension('/forecast/model/element') + '/forecast/model/element' + >>> f.nc_connectivity_dimension_groups() + () + + """ + return self._nc_clear_groups( + nc_get=self.nc_get_connectivity_dimension, + nc_set=self.nc_set_connectivity_dimension, + nc_groups=self.nc_connectivity_dimension_groups, + ) diff --git a/cfdm/mixin/propertiesdata.py b/cfdm/mixin/propertiesdata.py index d2d91fba2..7efc5d98b 100644 --- a/cfdm/mixin/propertiesdata.py +++ b/cfdm/mixin/propertiesdata.py @@ -70,7 +70,16 @@ def __getitem__(self, indices): data = self.get_data(None, _fill_value=False) if data is not None: - new.set_data(data[indices], copy=False) + original_shape = data.shape + data = data[indices] + + if 0 in data.shape: + raise IndexError( + f"Index {indices!r} selects no elements from " + f"data with shape {original_shape}" + ) + + new.set_data(data, copy=False) return new diff --git a/cfdm/mixin/topology.py b/cfdm/mixin/topology.py index 11db8a0d1..00f757f99 100644 --- a/cfdm/mixin/topology.py +++ b/cfdm/mixin/topology.py @@ -107,7 +107,10 @@ def _normalise_cell_ids(cls, data, start_index, remove_empty_columns): move_missing_values = True if move_missing_values: - # Move missing values to the end of each row + # Move missing values to the end of each row. + # + # Note: this might reorder the each row (excluding the + # first column). data[:, 1:].sort(axis=1, endwith=True) if remove_empty_columns: diff --git a/cfdm/read_write/netcdf/netcdfread.py b/cfdm/read_write/netcdf/netcdfread.py index b94a15cb6..10579c027 100644 --- a/cfdm/read_write/netcdf/netcdfread.py +++ b/cfdm/read_write/netcdf/netcdfread.py @@ -11,7 +11,6 @@ from numbers import Integral from os.path import isdir, isfile, join from typing import Any -from uuid import uuid4 import numpy as np @@ -80,8 +79,6 @@ class Mesh: # The netCDF dimension spanned by the cells for each # location. E.g. {'node': 'nNodes', 'edge': 'nEdges'} ncdim: dict = field(default_factory=dict) - # A unique identifier for the mesh. E.g. 'df10184d806ef1a10f5035e' - mesh_id: Any = None class NetCDFRead(IORead): @@ -331,7 +328,9 @@ def ugrid_cell_connectivity_types(self): """ return { + "edge_edge_connectivity": "node", "face_face_connectivity": "edge", + "volume_volume_connectivity": "face", } def ugrid_mesh_topology_attributes(self): @@ -346,6 +345,7 @@ def ugrid_mesh_topology_attributes(self): """ return ( + "cf_role", "topology_dimension", "node_coordinates", "edge_coordinates", @@ -358,6 +358,7 @@ def ugrid_mesh_topology_attributes(self): "face_face_connectivity", "face_edge_connectivity", "edge_node_connectivity", + "edge_edge_connectivity", "edge_face_connectivity", "volume_node_connectivity", "volume_edge_connectivity", @@ -1979,7 +1980,7 @@ def read( " read_vars['parsed_aggregated_data'] =\n" f" {g['parsed_aggregated_data']}" ) # prgama: no cover - + # ------------------------------------------------------------ # List variables # @@ -2159,7 +2160,11 @@ def read( self._ugrid_parse_location_index_set(attributes) if debug: - logger.debug(f" UGRID meshes:\n {g['mesh']}") + from pprint import pformat + + logger.debug( + f" UGRID meshes:\n {pformat(g['mesh'])}" + ) # pragma: no cover # ------------------------------------------------------------ # Identify and parse all quantization container variables @@ -2231,17 +2236,39 @@ def read( all_fields_or_domains[ncvar] = field_or_domain # ------------------------------------------------------------ - # Create domain constructs from UGRID mesh topology variables + # Create domain constructs from UGRID mesh topology variables, + # but only for mesh locations that haven't been accounted for + # by data or domain variables. # ------------------------------------------------------------ if domain and g["UGRID_version"] is not None: locations = ("node", "edge", "face") - for ncvar in g["variables"]: - if ncvar not in g["mesh"]: + for mesh_ncvar in g["variables"]: + if mesh_ncvar not in g["mesh"]: continue for location in locations: + # If any existing field or domain used this + # mesh/location combination, then we don't need to + # create a another new domain for it. + create_domain = True + for ncvar in all_fields_or_domains: + attributes = g["variable_attributes"].get(ncvar) + if attributes is None: + continue + + if ( + attributes.get("mesh") == mesh_ncvar + and attributes.get("location") == location + ): + create_domain = False + break + + if not create_domain: + continue + + # Still here? Create a new domain. mesh_domain = self._create_field_or_domain( - ncvar, + mesh_ncvar, domain=domain, location=location, ) @@ -2316,6 +2343,7 @@ def read( [ncvar for ncvar in sorted(g["do_not_create_field"])] ) ) # pragma: no cover + logger.info( " Unreferenced netCDF variables:\n " + "\n ".join(unreferenced_variables) @@ -3929,10 +3957,18 @@ def _create_field_or_domain( g["dataset_compliance"][field_ncvar]["dimensions"] = dimensions g["dataset_compliance"][field_ncvar].setdefault("non-compliance", {}) - logger.info( - " Converting netCDF variable " - f"{field_ncvar}({', '.join(dimensions)}) to a {construct_type}:" - ) # pragma: no cover + if mesh_topology: + logger.info( + " Converting netCDF mesh topology variable " + f"{field_ncvar}({', '.join(dimensions)}) to a " + f"{construct_type}:" + ) # pragma: no cover + else: + logger.info( + " Converting netCDF variable " + f"{field_ncvar}({', '.join(dimensions)}) to a " + f"{construct_type}:" + ) # pragma: no cover # ------------------------------------------------------------ # Combine the global and group properties with the data @@ -4737,12 +4773,13 @@ def _create_field_or_domain( # ------------------------------------------------------------ # Add a domain topology construct derived from UGRID + # (CF>=1.11) # ------------------------------------------------------------ if ugrid: domain_topology = mesh.domain_topologies.get(location) if domain_topology is not None: logger.detail( - " [m] Inserting " + f" [m] Inserting {location} " f"{domain_topology.__class__.__name__} with data shape " f"{self.implementation.get_data_shape(domain_topology)}" ) # pragma: no cover @@ -4759,13 +4796,14 @@ def _create_field_or_domain( # ------------------------------------------------------------ # Add a cell connectivity construct derived from UGRID + # (CF>=1.11) # ------------------------------------------------------------ if ugrid: for cell_connectivity in mesh.cell_connectivities.get( location, () ): logger.detail( - " [n] Inserting " + f" [n] Inserting {location} " f"{cell_connectivity.__class__.__name__} with data shape " f"{self.implementation.get_data_shape(cell_connectivity)}" ) # pragma: no cover @@ -4780,10 +4818,6 @@ def _create_field_or_domain( ncvar = self.implementation.nc_get_variable(cell_connectivity) ncvar_to_key[ncvar] = key - if ugrid: - # Set the mesh identifier - self.implementation.set_mesh_id(f, mesh.mesh_id) - # ------------------------------------------------------------ # Add coordinate reference constructs from formula_terms # properties @@ -9688,6 +9722,7 @@ def _ugrid_parse_mesh_topology(self, mesh_ncvar, attributes): "face_face_connectivity", "face_edge_connectivity", "edge_node_connectivity", + "edge_edge_connectivity", "edge_face_connectivity", "volume_node_connectivity", "volume_edge_connectivity", @@ -9703,7 +9738,6 @@ def _ugrid_parse_mesh_topology(self, mesh_ncvar, attributes): mesh = Mesh( mesh_ncvar=mesh_ncvar, mesh_attributes=attributes, - mesh_id=uuid4().hex, ) locations = ("node", "edge", "face") @@ -9847,7 +9881,6 @@ def _ugrid_parse_location_index_set(self, parent_attributes): location_index_set_attributes=location_index_set_attributes, location=location, index_set=index_set, - mesh_id=uuid4().hex, ) def _ugrid_create_auxiliary_coordinates( @@ -10139,9 +10172,7 @@ def _ugrid_create_domain_topology(self, parent_ncvar, f, mesh, location): # Create data if cell == "point": - properties["long_name"] = ( - "Maps every point to its connected points" - ) + properties["long_name"] = "Maps every node to its connected nodes" indices, kwargs = self._create_netcdfarray(connectivity_ncvar) n_nodes = self.read_vars["internal_dimension_sizes"][ mesh.ncdim[location] @@ -10190,6 +10221,12 @@ def _ugrid_create_domain_topology(self, parent_ncvar, f, mesh, location): # Apply a location index set domain_topology = domain_topology[index_set] + # Store the netCDF connectivity dimension name + connectivity_ncdim = self._ncdim_abspath( + self.read_vars["variable_dimensions"][connectivity_ncvar][-1] + ) + domain_topology.nc_set_connectivity_dimension(connectivity_ncdim) + return domain_topology def _ugrid_create_cell_connectivities( @@ -10224,7 +10261,7 @@ def _ugrid_create_cell_connectivities( location. May be an empty list. """ - if location != "face": + if location not in ("face", "edge"): return [] attributes = mesh.mesh_attributes diff --git a/cfdm/read_write/netcdf/netcdfwrite.py b/cfdm/read_write/netcdf/netcdfwrite.py index ec3b803d1..0314c0dee 100644 --- a/cfdm/read_write/netcdf/netcdfwrite.py +++ b/cfdm/read_write/netcdf/netcdfwrite.py @@ -21,6 +21,7 @@ ZARR_FMTS, ) from .netcdfread import NetCDFRead +from .netcdfwrite_ugrid import NetCDFWriteUgrid logger = logging.getLogger(__name__) @@ -35,8 +36,8 @@ class AggregationError(Exception): pass -class NetCDFWrite(IOWrite): - """A container for writing Fields to a netCDF dataset. +class NetCDFWrite(NetCDFWriteUgrid, IOWrite): + """A container for writing Fields to a dataset. NetCDF3, netCDF4 and Zarr output formats are supported. @@ -705,14 +706,14 @@ def _write_dimension( except RuntimeError as error: message = ( "Can't create unlimited dimension " - f"in {g['netcdf'].data_model} dataset ({error})." + f"in {g['dataset'].data_model} dataset ({error})." ) error = str(error) if error == "NetCDF: NC_UNLIMITED size already in use": raise RuntimeError( message - + f" In a {g['netcdf'].data_model} dataset only " + + f" In a {g['dataset'].data_model} dataset only " "one unlimited dimension is allowed. Consider " "using a netCDF4 format." ) @@ -726,7 +727,7 @@ def _write_dimension( except RuntimeError as error: raise RuntimeError( f"Can't create size {size} dimension {ncdim!r} in " - f"{g['netcdf'].data_model} dataset ({error})" + f"{g['dataset'].data_model} dataset ({error})" ) g["dimensions"].add(ncdim) @@ -1014,7 +1015,7 @@ def _write_scalar_data(self, f, value, ncvar): return ncvar def _create_geometry_container(self, field): - """Create a geometry container variable in the dataset. + """Create a geometry container. .. versionadded:: (cfdm) 1.8.0 @@ -1348,7 +1349,7 @@ def _write_bounds( # CF>=1.8 and we have geometry bounds, which are dealt # with separately # -------------------------------------------------------- - extra = self._write_node_coordinates( + extra = self._write_geometry_node_coordinates( f, coord, coord_ncvar, coord_ncdimensions ) return extra @@ -1467,10 +1468,10 @@ def _write_bounds( return extra - def _write_node_coordinates( + def _write_geometry_node_coordinates( self, f, coord, coord_ncvar, coord_ncdimensions ): - """Create a node coordinates dataset variable. + """Create a geometry node coordinates dataset variable. This will create: @@ -2263,7 +2264,6 @@ def _write_auxiliary_coordinate(self, f, key, coord, coordinates): ) else: ncvar = self._create_variable_name(coord, default="auxiliary") - # TODO: move setting of bounds ncvar to here - why? # If this auxiliary coordinate has bounds then create @@ -2281,11 +2281,16 @@ def _write_auxiliary_coordinate(self, f, key, coord, coordinates): self.implementation.get_data_axes(f, key), extra=extra, ) + else: + # There is no data (although there might be + # bounds, though), so there can't be an auxiliary + # coordinate variable. + ncvar = None - g["key_to_ncvar"][key] = ncvar g["key_to_ncdims"][key] = ncdimensions if ncvar is not None: + g["key_to_ncvar"][key] = ncvar coordinates.append(ncvar) return coordinates @@ -3746,14 +3751,13 @@ def _write_field_or_domain( ).get("grid_mapping_name", False) ] - # Check if the field or domain has a domain topology construct - # (CF>=1.11) - ugrid = self.implementation.has_domain_topology(f) - if ugrid: - raise NotImplementedError( - "Can't yet write UGRID datasets. " - "This feature is coming soon ..." - ) + # Initialise the dictionary of the field/domain's normalised + # domain topologies + g["normalised_domain_topologies"] = {} + + # Initialise the dictionary of the field/domain's normalised + # cell connectivities + g["normalised_cell_connectivities"] = {} field_coordinates = self.implementation.get_coordinates(f) @@ -4394,6 +4398,18 @@ def _write_field_or_domain( ).items() ] + # ------------------------------------------------------------ + # Domain topology variables (CF>=1.11) + # ------------------------------------------------------------ + for key, dt in self.implementation.get_domain_topologies(f).items(): + self._ugrid_write_domain_topology(f, key, dt) + + # ------------------------------------------------------------ + # Cell connectivity variables (CF>=1.11) + # ------------------------------------------------------------ + for key, cc in self.implementation.get_cell_connectivities(f).items(): + self._ugrid_write_cell_connectivity(f, key, cc) + # ------------------------------------------------------------ # Create the data/domain dataset variable # ------------------------------------------------------------ @@ -4402,7 +4418,7 @@ def _write_field_or_domain( else: default = "domain" - ncvar = self._create_variable_name(f, default=default) + field_ncvar = self._create_variable_name(f, default=default) ncdimensions = data_ncdimensions @@ -4413,7 +4429,7 @@ def _write_field_or_domain( cell_measures = " ".join(cell_measures) logger.debug( " Writing cell_measures attribute to " - f"variable {ncvar}: {cell_measures!r}" + f"variable {field_ncvar}: {cell_measures!r}" ) # pragma: no cover extra["cell_measures"] = cell_measures @@ -4423,7 +4439,7 @@ def _write_field_or_domain( coordinates = " ".join(coordinates) logger.info( " Writing coordinates attribute to " - f"variable {ncvar}: {coordinates!r}" + f"variable {field_ncvar}: {coordinates!r}" ) # pragma: no cover extra["coordinates"] = coordinates @@ -4433,7 +4449,7 @@ def _write_field_or_domain( grid_mapping = " ".join(grid_mapping) logger.info( " Writing grid_mapping attribute to " - f"variable {ncvar}: {grid_mapping!r}" + f"variable {field_ncvar}: {grid_mapping!r}" ) # pragma: no cover extra["grid_mapping"] = grid_mapping @@ -4444,7 +4460,7 @@ def _write_field_or_domain( ancillary_variables = re.sub(r"\s+", " ", ancillary_variables) logger.info( " Writing ancillary_variables attribute to " - f"variable {ncvar}: {ancillary_variables!r}" + f"variable {field_ncvar}: {ancillary_variables!r}" ) # pragma: no cover extra["ancillary_variables"] = ancillary_variables @@ -4481,7 +4497,7 @@ def _write_field_or_domain( cell_methods = " ".join(cell_methods_strings) logger.info( " Writing cell_methods attribute to " - f"variable {ncvar}: {cell_methods}" + f"variable {field_ncvar}: {cell_methods}" ) # pragma: no cover extra["cell_methods"] = cell_methods @@ -4497,6 +4513,24 @@ def _write_field_or_domain( ) extra["geometry"] = gc_ncvar + # ------------------------------------------------------------ + # UGRID mesh (CF>=1.11) + # + # Note this must be done *after* domain topologies and cell + # connectivities have been written to the dataset, as we'll + # need their netCDF variable names. + # ------------------------------------------------------------ + if g["output_version"] >= g["CF-1.11"]: + mesh_ncvar = self._ugrid_get_mesh_ncvar(f) + if mesh_ncvar is not None: + extra["mesh"] = mesh_ncvar + + location = f.domain_topology().get_cell() + if location == "point": + location = "node" + + extra["location"] = location + # ------------------------------------------------------------ # Create a new data/domain dataset variable # ------------------------------------------------------------ @@ -4519,7 +4553,7 @@ def _write_field_or_domain( # automatically changed to () within the # _write_netcdf_variable method. CF-1.9 self._write_netcdf_variable( - ncvar, + field_ncvar, ncdimensions, f, self.implementation.get_data_axes(f, None), @@ -4533,7 +4567,7 @@ def _write_field_or_domain( if add_to_seen: seen[id_f] = { "variable": org_f, - "ncvar": ncvar, + "ncvar": field_ncvar, "ncdims": ncdimensions, } @@ -4712,7 +4746,7 @@ def _get_group(self, parent, groups): :Parameters: - parent: `netCDF4.Dateset` or `netCDF4.Group` or `Zarr.Group` + parent: `netCDF4.Dataset` or `netCDF4.Group` or `Zarr.Group` The group in which to find or create new group. groups: sequence of `str` @@ -5474,6 +5508,10 @@ def write( # by their output dataset variable names. # -------------------------------------------------------- "quantization": {}, + # -------------------------------------------------------- + # UGRID: + # -------------------------------------------------------- + "meshes": {}, } if mode not in ("w", "a", "r+"): @@ -5556,10 +5594,8 @@ def write( if mode == "a": if fmt == "ZARR3": - raise ValueError( - "Can't write with mode 'a' to a Zarr dataset" - ) - + raise ValueError("Can't write with mode 'a' to a Zarr dataset") + # First read in the fields from the existing dataset: effective_fields = self._NetCDFRead(self.implementation).read( dataset_name, netcdf_backend="netCDF4" @@ -5651,7 +5687,8 @@ def write( if mode == "w": # only one iteration required in this simple case return - elif mode == "a": # need another iteration to append after reading + + if mode == "a": # need another iteration to append after reading self.write_vars["dry_run"] = False self.write_vars["post_dry_run"] = True # i.e. follows a dry run @@ -5912,6 +5949,12 @@ def _file_io_iteration( for f in fields: self._write_field_or_domain(f) + # ------------------------------------------------------------ + # Now write any UGRID meshes referenced by any of the + # fields/domains (CF>=1.11) + # ------------------------------------------------------------ + self._ugrid_write_mesh_variables() + # ------------------------------------------------------------ # Write all of the buffered data to disk # ------------------------------------------------------------ diff --git a/cfdm/read_write/netcdf/netcdfwrite_ugrid.py b/cfdm/read_write/netcdf/netcdfwrite_ugrid.py new file mode 100644 index 000000000..e0bea4622 --- /dev/null +++ b/cfdm/read_write/netcdf/netcdfwrite_ugrid.py @@ -0,0 +1,1050 @@ +import logging + +import numpy as np + +logger = logging.getLogger(__name__) + + +class NetCDFWriteUgrid: + """Mixin class for writing UGRID meshes to a dataset. + + .. versionadded: (cfdm) NEXTVERSION + + """ + + def _ugrid_write_domain_topology(self, parent, key, domain_topology): + """Write a domain topology to a *_node_connectivity variable. + + If an equal domain topology has already been written to the + dataset then it is not re-written. + + .. versionadded: (cfdm) NEXTVERSION + + :Parameters: + + parent: `Field` or `Domain` or `None` + The parent Field or Domain. Set to `None` if there is + no parent (as could occur when for an "edge" + `domain_topology` that is derived from point-cell + domain topology). + + key: `str` or `None` + The internal identifier of the domain topology + construct. Set to `None` if *parent* is `None`. + + domain_topology: `DomainTopology` + The Domain Topology construct to be written. + + :Returns: + + `str` or `None` + The dataset variable name of the domain topology + object, or `None` if one wasn't written. + + **Examples** + + >>> ncvar = n._ugrid_write_domain_topology(f, 'domaintopology0', dt) + + """ + from cfdm.functions import integer_dtype + + g = self.write_vars + + # Normalise the array, so that its N node ids are 0, ..., N-1 + + domain_topology.normalise(inplace=True) + if key is not None: + g["normalised_domain_topologies"][key] = domain_topology + + cell = domain_topology.get_cell(None) + if cell == "point": + # There's no corresponding UGRID variable for "point" + # cells, so just return. + return None + + if cell == "volume": + raise NotImplementedError( + f"Can't write a UGRID mesh of volume cells for {parent!r}" + ) + + if cell not in ("face", "edge"): + raise ValueError( + f"{parent!r} has unknown domain topology cell type: " + f"{domain_topology!r}" + ) + + # Get the netCDF dimensions + size0, size1 = domain_topology.data.shape + if parent is not None: + # Get the number-of-cells dimension name from the parent + cells_ncdim = self._dataset_dimensions( + parent, key, domain_topology + ) + cells_ncdim = cells_ncdim[0] + else: + # Get the number-of-cells dimension name without reference + # to a parent + cells_ncdim = self._name( + f"{cell}s", dimsize=size0, role=f"ugrid_{cell}" + ) + + # Get the connectivity dimension name from the DomainTopology + connectivity_ncdim = domain_topology.nc_get_connectivity_dimension( + f"connectivity{size1}" + ) + connectivity_ncdim = self._name( + connectivity_ncdim, dimsize=size1, role="ugrid_connectivity" + ) + + ncdimensions = (cells_ncdim, connectivity_ncdim) + + if self._already_in_file(domain_topology, ncdimensions): + # This domain topology variable has been previously + # created, so no need to do so again. + ncvar = g["seen"][id(domain_topology)]["ncvar"] + else: + # This domain topology variable has not been previously + # created, so create it now. + + if cells_ncdim not in g["ncdim_to_size"]: + # Create a new number-of-cells netCDF dimension + self._write_dimension(cells_ncdim, parent, size=size0) + + if connectivity_ncdim not in g["ncdim_to_size"]: + # Create a new connectivity netCDF dimension + self._write_dimension(connectivity_ncdim, parent, size=size1) + + ncvar = self._create_variable_name( + domain_topology, default=f"{cell}_node_connectivity" + ) + + # Write as 32-bit integers if possible + dtype = integer_dtype(domain_topology.data.size) + domain_topology.data.dtype = dtype + + if parent is not None: + # Get domain axis keys from the parent + domain_axes = self.implementation.get_data_axes(parent, key) + else: + domain_axes = None + + self._write_netcdf_variable( + ncvar, + ncdimensions, + domain_topology, + domain_axes, + ) + + g["key_to_ncvar"][key] = ncvar + g["key_to_ncdims"][key] = ncdimensions + + return ncvar + + def _ugrid_write_cell_connectivity(self, f, key, cell_connectivity): + """Write a cell connectivity to a *_*_connectivity variable. + + If an equal cell connectivity has already been written to the + dataset then it is not re-written. + + .. versionadded: (cfdm) NEXTVERSION + + :Parameters: + + f: `Field` or `Domain` + The parent Field or Domain. + + key: `str` + The internal identifier of the cell connectivity + construct. + + cell_connectivity: `CellConnectivity` + The Cell Connectivity construct to be written. + + :Returns: + + `str` + The dataset variable name of the cell connectivity + object. + + **Examples** + + >>> ncvar = n._ugrid_write_cell_connectivity(f, 'cellconn0', cc) + + """ + from cfdm.functions import integer_dtype + + g = self.write_vars + + # Normalise the array, so that its N cell ids are 0, ..., N-1 + cell_connectivity.normalise(inplace=True) + g["normalised_cell_connectivities"][key] = cell_connectivity + + # Remove the first column, which (now that the array has been + # normalised) just contains the index of each row (0, ..., + # N-1) and is not stored in the dataset. + cell_connectivity = cell_connectivity[:, 1:] + + # Get the netCDF dimensions + size = cell_connectivity.data.shape[-1] + ncdimensions = self._dataset_dimensions(f, key, cell_connectivity) + connectivity_ncdim = cell_connectivity.nc_get_connectivity_dimension( + f"connectivity{size}" + ) + connectivity_ncdim = self._name( + connectivity_ncdim, dimsize=size, role="ugrid_connectivity" + ) + ncdimensions = ncdimensions + (connectivity_ncdim,) + + if self._already_in_file(cell_connectivity, ncdimensions): + # This cell_connectivity variable has been previously + # created, so no need to do so again. + ncvar = g["seen"][id(cell_connectivity)]["ncvar"] + else: + # This cell_connectivity variable has not been previously + # created, so create it now. + if connectivity_ncdim not in g["ncdim_to_size"]: + self._write_dimension(connectivity_ncdim, f, size=size) + + match cell_connectivity.get_connectivity("cell"): + case "edge": + cell = "face" + case "node": + cell = "edge" + case "face": + cell = "volume" + + ncvar = self._create_variable_name( + cell_connectivity, default=f"{cell}_{cell}_connectivity" + ) + + # Write as 32-bit integers if possible + dtype = integer_dtype(cell_connectivity.shape[0]) + cell_connectivity.data.dtype = dtype + + self._write_netcdf_variable( + ncvar, + ncdimensions, + cell_connectivity, + self.implementation.get_data_axes(f, key), + ) + + g["key_to_ncvar"][key] = ncvar + g["key_to_ncdims"][key] = ncdimensions + + return ncvar + + def _ugrid_write_node_coordinate(self, node_coord, ncdimensions): + """Write UGRID node coordinates to the dataset. + + .. versionadded:: (cfdm) NEXTVERSION + + :Parameters: + + node_coord: `AuxiliaryCoordinate` + The node coordinates to be written. + + ncdimensions: sequence of `str` + The netCDF dimension name of the node coordinates + dimension, e.g. ``('n_nodes',)``. + + :Returns: + + `str` + The netCDF variable name of the dateset node + coordinates variable. + + """ + g = self.write_vars + + already_in_file = self._already_in_file(node_coord, ncdimensions) + + if already_in_file: + ncvar = g["seen"][id(node_coord)]["ncvar"] + else: + ncvar = self._create_variable_name( + node_coord, default="node_coordinates" + ) + + # Create a new UGRID node coordinate variable + if self.implementation.get_data(node_coord, None) is not None: + self._write_netcdf_variable( + ncvar, + ncdimensions, + node_coord, + None, + ) + + return ncvar + + def _ugrid_get_mesh_ncvar(self, parent): + """Get the name of the netCDF mesh variable. + + .. versionadded:: (cfdm) NEXTVERSION + + :Parameters: + + parent: `Field` or `Domain` + The parent Field or Domain from which to get the mesh + description. + + :Returns: + + `str` or `None` + The name of the netCDF mesh variable that will store + the parent's UGRID mesh, or `None` if the parent is + not UGRID. + + """ + g = self.write_vars + + # Create a new mesh description for this parent + ncvar_new, mesh_new = self._ugrid_create_mesh(parent) + + if ncvar_new is None: + # Parent is not UGRID + return + + for ncvar, mesh in g["meshes"].items(): + if self._ugrid_linked_meshes(mesh, mesh_new): + # The mesh is either A) identical to another parent's + # mesh, or B) represents a different location (node, + # edge, face, or volume) of another parent's mesh. + # + # In both cases we can assign that other parent's mesh + # to the current parent; but in case B), we first + # update the other parent's mesh to include the new + # location. In case A), `_ugrid_update_mesh` makes no + # change to the other parent's mesh. + self._ugrid_update_mesh(mesh, mesh_new) + return ncvar + + # Still here? Then this parent's UGRID mesh is not the same + # as, nor linked to, any other parent's mesh, so we save save + # it as a new mesh. + g["meshes"][ncvar_new] = mesh_new + + return ncvar_new + + def _ugrid_create_mesh(self, parent): + """Create a mesh description from a parent Field or Domain. + + The mesh description is a dictionary that always has the + keys:: + + 'attributes' + 'topology_dimension' + 'node_coordinates' + 'sorted_edges' + + as well as a subset of the optional keys:: + + 'edge_coordinates' + 'face_coordinates' + 'volume_coordinates' + 'edge_node_connectivity' + 'face_node_connectivity' + 'volume_node_connectivity' + 'edge_edge_connectivity' + 'face_face_connectivity' + 'volume_volume_connectivity' + + After the mesh description is a dictionary has been created + with this method: + + * More of the optional keys might get added by + `_ugrid_update_mesh`. + * The 'attributes' dictionary might get updated by + `_ugrid_update_mesh`. + * The 'sorted_edges' dictionary might get updated by various + methods. + + For instance, consider the mesh description for the UGRID mesh + topology of face cells taken from ``cfdm.example_field(8)``. + In this case the 'node_coordinates' Auxiliary Coordinates are + derived from the face cell bounds:: + + {'attributes': + {'face_coordinates': ['Mesh2_face_x', 'Mesh2_face_y'], + 'face_face_connectivity': ['Mesh2_face_links'], + 'face_node_connectivity': ['Mesh2_face_nodes']}, + 'face_coordinates': + [, + ], + 'face_face_connectivity': + [], + 'face_node_connectivity': + [], + 'node_coordinates': + [, + ], + 'sorted_edges': + {}, + 'topology_dimension': + 2 + } + + For instance, consider the mesh description for the UGRID mesh + topology of edge cells taken from ``cfdm.example_field(9). In + this case the 'node_coordinates' Auxiliary Coordinates are + derived from the edge cell bounds:: + + {'attributes': + {'edge_coordinates': ['Mesh2_edge_x', 'Mesh2_edge_y'], + 'edge_edge_connectivity': ['Mesh2_edge_links'], + 'edge_node_connectivity': ['Mesh2_edge_nodes']}, + 'edge_coordinates': + [, + ], + 'edge_edge_connectivity': + [], + 'edge_node_connectivity': + [], + 'node_coordinates': + [, + ], + 'sorted_edges': + {}, + 'topology_dimension': + 1 + } + + For instance, consider the mesh description for the UGRID mesh + topology of node cells taken from ``cfdm.example_field(10). In + this case the 'node_coordinates' Auxiliary Coordinates are + explicitly defined by the point cell locations:: + + {'attributes': + {'edge_node_connectivity': [], + 'node_coordinates': ['Mesh2_node_x', 'Mesh2_node_y'], + 'node_node_connectivity': []}, + 'node_coordinates': + [, + ], + 'node_node_connectivity': + [], + 'sorted_edges': + {}, + 'topology_dimension': + 0 + } + + .. versionadded:: (cfdm) NEXTVERSION + + :Parameters: + + parent: `Field` or `Domain` + The parent Field or Domain from which to create the + mesh description. + + :Returns: + + (`str`, `dict`) or (`None`, `None`) + The name of the netCDF mesh variable that will store + the parent's UGRID mesh, and the mesh description + itself. If the parent is not UGRID, then both are + returned as `None`. + + """ + g = self.write_vars + + # Get the dictionary of normalised domain topology constructs + domain_topologies = g["normalised_domain_topologies"] + if not domain_topologies: + # Not UGRID + return None, None + + # Initialise the mesh description. + # + # This always includes the sub-dictionary 'attributes', which + # contains the netCDF names of mesh-related variables; and the + # sub-dictionary 'sorted_edges', which is a cache of the + # sorted unique edges implied by the domain topologies. + mesh = {"attributes": {}, "sorted_edges": {}} + + if len(domain_topologies) > 1: + raise ValueError( + f"Can't write a UGRID mesh for {parent!r} that has multiple " + "domain topology constructs: " + f"{tuple(domain_topologies.values())}" + ) + + # Get the UGRID domain axis + key, domain_topology = domain_topologies.popitem() + + ugrid_axis = self.implementation.get_data_axes(parent, key)[0] + + # Get the dataset variable name of the domain topology + # construct + ncvar_cell_node_connectivity = g["key_to_ncvar"].get(key, []) + if ncvar_cell_node_connectivity != []: + ncvar_cell_node_connectivity = [ncvar_cell_node_connectivity] + + # Get the 1-d auxiliary coordinates that span the UGRID axis, + # and their dataset variable names + cell_coordinates = self.implementation.get_auxiliary_coordinates( + parent, axes=(ugrid_axis,), exact=True + ) + ncvar_cell_coordinates = [ + g["key_to_ncvar"][key] + for key in cell_coordinates + if key in g["key_to_ncvar"] + ] + + # ------------------------------------------------------------ + # Populate the output mesh description + # ------------------------------------------------------------ + cell = domain_topology.get_cell(None) + if cell == "point": + mesh.update( + { + "topology_dimension": 0, + "node_coordinates": list(cell_coordinates.values()), + "node_node_connectivity": [domain_topology], + } + ) + mesh["attributes"].update( + { + "node_coordinates": ncvar_cell_coordinates, + "edge_node_connectivity": ncvar_cell_node_connectivity, + "node_node_connectivity": [], + } + ) + else: + match cell: + case "face": + topology_dimension = 2 + case "edge": + topology_dimension = 1 + case "volume": + topology_dimension = 3 + raise NotImplementedError( + "Can't write a UGRID mesh of volume cells for " + f"{parent!r}" + ) + case _: + raise ValueError( + f"{parent!r} has unknown domain topology cell type: " + f"{domain_topology!r}" + ) + + node_coordinates = [] + index = None + for key, c in cell_coordinates.items(): + bounds = c.get_bounds(None) + if bounds is None or bounds.get_data(None) is None: + raise ValueError( + f"Can't write a UGRID mesh of {cell!r} cells for " + f"{parent!r} when {c!r} has no coordinate bounds data" + ) + + if index is None: + # Create the array index that will extract, in the + # correct order, the node coordinates from the + # flattened cell bounds, i.e. so that the first + # node coordinate has node id 0, the second has + # node id 1, etc. + node_ids, index = np.unique( + domain_topology, return_index=True + ) + if node_ids[-1] is np.ma.masked: + # Remove the element that corresponds to + # missing data (which is always at the end of + # the `np.unique` outputs) + index = index[:-1] + + # Create, from the cell bounds, an Auxiliary + # Coordinate that contains the unique node + # coordinates. + coords = self.implementation.initialise_AuxiliaryCoordinate( + data=bounds.data.flatten()[index], + properties=c.properties(), + ) + + # Persist the node coordinates into memory, because + # it's possible that we'll need to compare them more + # than once with the node coordinates of other mesh + # descriptions (in `_ugrid_linked_meshes`), in + # addition to writing the coordinate themselves to + # disk. + coords.persist(inplace=True) + + node_coordinates.append(coords) + + del index + + mesh.update( + { + "topology_dimension": topology_dimension, + "node_coordinates": node_coordinates, + f"{cell}_node_connectivity": [domain_topology], + } + ) + mesh["attributes"][ + f"{cell}_node_connectivity" + ] = ncvar_cell_node_connectivity + + if ncvar_cell_coordinates: + key = f"{cell}_coordinates" + mesh[key] = cell_coordinates + mesh["attributes"][key] = ncvar_cell_coordinates + + # Add mesh description keys for normalised cell connectivities + for cc_key, cell_connectivity in g[ + "normalised_cell_connectivities" + ].items(): + connectivity = cell_connectivity.get_connectivity(None) + if not ( + (connectivity, cell) == ("edge", "face") + or (connectivity, cell) == ("node", "edge") + or (connectivity, cell) == ("face", "volume") + ): + raise ValueError( + f"{parent!r} has invalid UGRID cell connectivity type " + f"for {cell!r} cells: {cell_connectivity!r}" + ) + + key = f"{cell}_{cell}_connectivity" + mesh[key] = [cell_connectivity] + mesh["attributes"][key] = [g["key_to_ncvar"][cc_key]] + + # Get the dataset mesh variable name + ncvar = parent.nc_get_mesh_variable("mesh") + ncvar = self._name(ncvar) + + return ncvar, mesh + + def _ugrid_linked_meshes(self, mesh, mesh1): + """Ascertain if two meshes are linked. + + Meshes are linked if they represent different locations of the + same UGRID mesh. + + .. versionadded:: (cfdm) NEXTVERSION + + :Parameters: + + mesh: `dict` + The two mesh dictionaries to be compared. + + mesh1: `dict` + The two mesh dictionaries to be compared. + + :Returns: + + `bool` + + """ + # Find the relevant keys that are common to both meshes + keys = ( + "node_coordinates", + "edge_coordinates", + "face_coordinates", + "volume_coordinates", + "node_node_connectivity", + "edge_node_connectivity", + "face_node_connectivity", + "volume_node_connectivity", + "edge_edge_connectivity", + "face_face_connectivity", + "volume_volume_connectivity", + ) + common_keys = [k for k in keys if k in mesh and k in mesh1] + + # Check the common keys for equality + # + # A necessary condition for meshes being linked is that these + # common keys have identical values. + for key in common_keys: + if len(mesh[key]) != len(mesh1[key]): + # Different numbers of constructs for the same key => + # the meshes are not linked + return False + + mesh1_key = mesh1[key][:] + for c in mesh[key]: + found_match = False + for i, c1 in enumerate(mesh1_key): + if c.equals(c1): + # Matching construct pair + found_match = True + mesh1_key.pop(i) + break + + if not found_match: + # A 'mesh1' construct is different to all those in + # 'mesh' => the meshes are not linked + return False + + # Still here? Then all of the keys common to both meshes have + # equal values. + + # Now check the non-common connectivity keys for consistency + # (as opposed to equality) + locations = ("edge", "node", "face", "volume") + location_mesh = {} + for location in locations: + key = f"{location}_node_connectivity" + if key in common_keys: + continue + + if key in mesh1: + location_mesh[location] = mesh1 + break + + for location in locations: + key = f"{location}_node_connectivity" + if key in common_keys or key in location_mesh: + continue + + if key in mesh: + location_mesh[location] = mesh + break + + if len(location_mesh) == 2: + # 'location_mesh' has two keys, one for each mesh, and + # each key represents a domain topology that is not + # present in the other mesh. + # + # Each pair of domain topology cell types needs special + # treatment. + if set(location_mesh) == set(("edge", "face")): + if not self._ugrid_check_edge_face(**location_mesh): + return False + + elif set(location_mesh) == set(("node", "edge")): + if not self._ugrid_check_node_edge(**location_mesh): + return False + + elif set(location_mesh) == set(("node", "face")): + if not self._ugrid_check_node_face(**location_mesh): + return False + + elif "volume" in location_mesh: + raise NotImplementedError( + "Can't write a UGRID mesh of volume cells" + ) + + # Still here? Then 'mesh' and 'mesh1' are parts of the same + # uber-mesh. + return True + + def _ugrid_check_node_edge(self, node, edge): + """Whether or not the nodes imply the edges, and vice versa. + + .. versionadded:: (cfdm) NEXTVERSION + + :Parameters: + + node: `dict` + The mesh dictionary of the nodes. + + edge: `dict` + The mesh dictionary of the edges. + + :Returns: + + `bool` + The result. + + """ + # Find the set of unique edges that are implied by the nodes + node_edges = node["sorted_edges"].get("node_node_connectivity") + if node_edges is None: + node_edges = node["node_node_connectivity"][0].to_edge(sort=True) + node["sorted_edges"]["node_node_connectivity"] = node_edges + node["sorted_edges"]["edge_node_connectivity"] = node_edges + + edges = edge["sorted_edges"].get("edge_node_connectivity") + if edges is None: + edges = edge["edge_node_connectivity"][0].sort() + edge["sorted_edges"]["edge_node_connectivity"] = edges + + # Return True if the unique edges of the faces are identical + # to the given edges + if node_edges.data.shape != edges.data.shape: + return False + + return bool((node_edges.data == edges.data).all()) + + def _ugrid_check_edge_face(self, edge, face): + """Whether or not the edges imply the faces, and vice versa. + + .. versionadded:: (cfdm) NEXTVERSION + + :Parameters: + + edge: `dict` + The mesh dictionary of the edges. + + face: `dict` + The mesh dictionary of the faces. + + :Returns: + + `bool` + The result. + + """ + # Fast checks that are sufficient (but not necessary) + # conditions for the edges and faces being incompatible + edges = edge["edge_node_connectivity"][0] + faces = face["face_node_connectivity"][0] + if edges.size > faces.data.size: + # There are more edges than the maximum that could be + # defined by the faces + return False + + if edges.size < faces.size + 2: + # There are fewer edges than the minimum that could be + # defined by the faces + return False + + # Still here? + edges = edge["sorted_edges"].get("edge_node_connectivity") + if edges is None: + edges = edge["edge_node_connectivity"][0].sort() + edge["sorted_edges"]["edge_node_connectivity"] = edges + + # Find the set of unique edges that are implied by the faces + face_edges = face["sorted_edges"].get("face_node_connectivity") + if face_edges is None: + n_nodes = face["node_coordinates"][0].size + face_edges = face["face_node_connectivity"][0].to_edge( + sort=True, face_nodes=range(n_nodes) + ) + face["sorted_edges"]["face_node_connectivity"] = face_edges + face["sorted_edges"]["edge_node_connectivity"] = face_edges + + # Return True if the unique edges of the faces are identical + # to the given edges + if face_edges.data.shape != edges.data.shape: + return False + + return bool((face_edges.data == edges.data).all()) + + def _ugrid_check_node_face(self, node, face): + """Whether or not the nodes imply the faces, and vice versa. + + .. versionadded:: (cfdm) NEXTVERSION + + :Parameters: + + node: `dict` + The mesh dictionary of the nodes. + + face: `dict` + The mesh dictionary of the faces. + + :Returns: + + `bool` + The result. + + """ + # Find the set of unique edges that are implied by the nodes + node_edges = node["sorted_edges"].get("node_node_connectivity") + if node_edges is None: + node_edges = node["node_node_connectivity"][0].to_edge(sort=True) + node["sorted_edges"]["node_node_connectivity"] = node_edges + node["sorted_edges"]["edge_node_connectivity"] = node_edges + + # Find the set of unique edges that are implied by the faces + face_edges = face["sorted_edges"].get("face_node_connectivity") + if face_edges is None: + n_nodes = face["node_coordinates"][0].size + face_edges = face["face_node_connectivity"][0].to_edge( + face_nodes=range(n_nodes), sort=True + ) + face["sorted_edges"]["face_node_connectivity"] = face_edges + face["sorted_edges"]["edge_node_connectivity"] = face_edges + + # Return True if the unique edges of the faces are identical + # to the given edges + if face_edges.data.shape != node_edges.data.shape: + return False + + return bool((face_edges.data == node_edges.data).all()) + + def _ugrid_update_mesh(self, mesh, mesh1): + """Update an original mesh with a linked mesh. + + Elements of the linked mesh which are not in the original mesh + are copied to the original mesh. + + The `_ugrid_linked_meshes` method is used to ascertain if two + meshes are linked. + + .. versionadded:: (cfdm) NEXTVERSION + + :Parameters: + + mesh: `dict` + The original mesh to be updated in-place. If `mesh` + and `mesh1` are identical then `mesh` is unchanged. + + mesh1: `dict` + The linked mesh to update from. + + :Returns: + + `None` + + """ + # Update the topology dimension + mesh["topology_dimension"] = max( + mesh["topology_dimension"], mesh1["topology_dimension"] + ) + + for key, value in mesh1.items(): + if key not in mesh: + # This 'mesh1' key is not in 'mesh', so copy it from + # 'mesh1'. + # + # Note: any such key will always have a `list` value. + mesh[key] = value.copy() + mesh["attributes"][key] = mesh1["attributes"][key].copy() + + for key, value in mesh1["sorted_edges"].items(): + if key not in mesh["sorted_edges"]: + # This mesh1["sorted_edges"] key is not in + # mesh["sorted_edges"], so copy it from + # mesh1["sorted_edges"]. + # + # Note: any such key will always have a + # `DomainTopology` value. + mesh["sorted_edges"][key] = value.copy() + + # If applicable, make sure that the node coordinates and their + # netCDF variable names are defined by the point-cell domain, + # rather than being inferred by one of the + # edge/face/volume-cell domains. + if "node_node_connectivity" in mesh1: + key = "node_coordinates" + mesh[key] = mesh1[key].copy() + mesh["attributes"][key] = mesh1["attributes"][key].copy() + + def _ugrid_write_mesh_variables(self): + """Write mesh variables to the dataset. + + This is done after all Fields and Domains have been written to + the dataset. + + All UGRID CF-netCDF data and domain variables in the dataset + will already have the correct mesh variable name stored in + their 'mesh' attributes. + + The mesh variables are defined by `self.write_vars['meshes']`. + + .. versionadded:: (cfdm) NEXTVERSION + + :Returns: + + `None` + + """ + g = self.write_vars + + for mesh_ncvar, mesh in g["meshes"].items(): + # -------------------------------------------------------- + # Create the mesh variable attributes. + # + # E.g. the 'mesh' dictionary + # + # {'attributes': + # {'face_coordinates': ['Mesh2_face_x', 'Mesh2_face_y'], + # 'face_face_connectivity': ['Mesh2_face_links'], + # 'face_node_connectivity': ['Mesh2_face_nodes']}, + # 'topology_dimension': + # 2, + # 'face_coordinates': + # [, + # ], + # 'face_node_connectivity': + # [], + # 'node_coordinates': + # [, + # ]} + # + # could give mesh variable attributes: + # + # {'topology_dimension': 2, + # 'node_coordinates': 'longitude latitude', + # 'face_coordinates': 'Mesh2_face_x Mesh2_face_y', + # 'face_node_connectivity': 'Mesh2_face_nodes'} + # -------------------------------------------------------- + attributes = { + "cf_role": "mesh_topology", + "topology_dimension": mesh["topology_dimension"], + } + + # Convert non-empty lists of constructs to a + # space-separated dataset variable name string + # + # E.g. [] -> 'Mesh2_face_links' + # E.g. [, ] -> 'x y' + for key, value in mesh["attributes"].items(): + if value: + attributes[key] = " ".join(value) + + # If the dataset variable names for the node coordinates + # (which have to exist for all meshes) have not been + # defined, then it's because the node coordinates have not + # yet been written to the dataset (which in turn is + # because there were no node-location field or domain + # constructs being written). So, let's write them now, and + # get their variable names. + if "node_coordinates" not in attributes: + # Create a new node dimension in the same group as the + # mesh variable + ncdim = "nodes" + mesh_group = self._groups(mesh_ncvar) + if mesh_group: + ncdim = mesh_group + ncdim + + n_nodes = mesh["node_coordinates"][0].size + ncdim = self._name(ncdim, dimsize=n_nodes, role="nodes") + if ncdim not in g["ncdim_to_size"]: + self._write_dimension(ncdim, None, size=n_nodes) + + # Write the node coordinates to the dataset + ncvars = [ + self._ugrid_write_node_coordinate(nc, (ncdim,)) + for nc in mesh["node_coordinates"] + ] + attributes["node_coordinates"] = " ".join(ncvars) + + # For a point-cell domain topology + # (i.e. 'topology_dimension' is currently 0), write the + # implied edge_node_connectivity variable to the dataset. + if not mesh["topology_dimension"]: + edges = mesh["sorted_edges"].get("node_node_connectivity") + if edges is None: + edges = mesh["node_node_connectivity"][0].to_edge( + sort=True + ) + + ncvar = self._ugrid_write_domain_topology(None, None, edges) + if ncvar is not None: + attributes["edge_node_connectivity"] = ncvar + + # Set topology dimension to 1, now that we've + # included an edge_node_connectivity. + attributes["topology_dimension"] = 1 + + # -------------------------------------------------------- + # Create the mesh variable and set its attributes + # -------------------------------------------------------- + logger.debug( + f" Writing UGRID mesh variable: {mesh_ncvar}\n" + f" {mesh}" + ) # pragma: no cover + + kwargs = { + "varname": mesh_ncvar, + "datatype": "S1", + "endian": g["endian"], + } + kwargs.update(g["netcdf_compression"]) + + self._createVariable(**kwargs) + self._set_attributes(attributes, mesh_ncvar) diff --git a/cfdm/test/create_test_files.py b/cfdm/test/create_test_files.py index 6f1e0db66..a78cb39c1 100644 --- a/cfdm/test/create_test_files.py +++ b/cfdm/test/create_test_files.py @@ -2229,6 +2229,150 @@ def _make_ugrid_2(filename): return filename +def _make_ugrid_3(filename): + """Create a UGRID mesh topology and no fields/domains.""" + n = netCDF4.Dataset(filename, "w") + + n.Conventions = f"CF-{VN}" + + n.createDimension("nMesh3_node", 7) + n.createDimension("nMesh3_edge", 9) + n.createDimension("nMesh3_face", 3) + n.createDimension("connectivity2", 2) + n.createDimension("connectivity4", 4) + n.createDimension("connectivity5", 5) + + Mesh3 = n.createVariable("Mesh3", "i4", ()) + Mesh3.cf_role = "mesh_topology" + Mesh3.topology_dimension = 2 + Mesh3.node_coordinates = "Mesh3_node_x Mesh3_node_y" + Mesh3.face_node_connectivity = "Mesh3_face_nodes" + Mesh3.edge_node_connectivity = "Mesh3_edge_nodes" + Mesh3.face_dimension = "nMesh3_face" + Mesh3.edge_dimension = "nMesh3_edge" + Mesh3.face_face_connectivity = "Mesh3_face_links" + Mesh3.edge_edge_connectivity = "Mesh3_edge_links" + + # Node + Mesh3_node_x = n.createVariable("Mesh3_node_x", "f4", ("nMesh3_node",)) + Mesh3_node_x.standard_name = "longitude" + Mesh3_node_x.units = "degrees_east" + Mesh3_node_x[...] = [-45, -43, -45, -43, -45, -43, -40] + + Mesh3_node_y = n.createVariable("Mesh3_node_y", "f4", ("nMesh3_node",)) + Mesh3_node_y.standard_name = "latitude" + Mesh3_node_y.units = "degrees_north" + Mesh3_node_y[...] = [35, 35, 33, 33, 31, 31, 34] + + Mesh3_edge_nodes = n.createVariable( + "Mesh3_edge_nodes", "i4", ("nMesh3_edge", "connectivity2") + ) + Mesh3_edge_nodes.long_name = "Maps every edge to its two nodes" + Mesh3_edge_nodes[...] = [ + [1, 6], + [3, 6], + [3, 1], + [0, 1], + [2, 0], + [2, 3], + [2, 4], + [5, 4], + [3, 5], + ] + + # Face + Mesh3_face_x = n.createVariable( + "Mesh3_face_x", "f8", ("nMesh3_face",), fill_value=-99 + ) + Mesh3_face_x.standard_name = "longitude" + Mesh3_face_x.units = "degrees_east" + Mesh3_face_x[...] = [-44, -44, -42] + + Mesh3_face_y = n.createVariable( + "Mesh3_face_y", "f8", ("nMesh3_face",), fill_value=-99 + ) + Mesh3_face_y.standard_name = "latitude" + Mesh3_face_y.units = "degrees_north" + Mesh3_face_y[...] = [34, 32, 34] + + Mesh3_face_nodes = n.createVariable( + "Mesh3_face_nodes", + "i4", + ("nMesh3_face", "connectivity4"), + fill_value=-99, + ) + Mesh3_face_nodes.long_name = "Maps every face to its corner nodes" + Mesh3_face_nodes[...] = [[2, 3, 1, 0], [4, 5, 3, 2], [6, 1, 3, -99]] + + Mesh3_face_links = n.createVariable( + "Mesh3_face_links", + "i4", + ("nMesh3_face", "connectivity4"), + fill_value=-99, + ) + Mesh3_face_links.long_name = "neighbour faces for faces" + Mesh3_face_links[...] = [ + [1, 2, -99, -99], + [0, -99, -99, -99], + [0, -99, -99, -99], + ] + + # Edge + Mesh3_edge_x = n.createVariable( + "Mesh3_edge_x", "f8", ("nMesh3_edge",), fill_value=-99 + ) + Mesh3_edge_x.standard_name = "longitude" + Mesh3_edge_x.units = "degrees_east" + Mesh3_edge_x[...] = [-41.5, -41.5, -43, -44, -45, -44, -45, -44, -43] + + Mesh3_edge_y = n.createVariable( + "Mesh3_edge_y", "f8", ("nMesh3_edge",), fill_value=-99 + ) + Mesh3_edge_y.standard_name = "latitude" + Mesh3_edge_y.units = "degrees_north" + Mesh3_edge_y[...] = [34.5, 33.5, 34, 35, 34, 33, 32, 31, 32] + + Mesh3_edge_links = n.createVariable( + "Mesh3_edge_links", + "i4", + ("nMesh3_edge", "connectivity5"), + fill_value=-99, + ) + Mesh3_edge_links.long_name = "neighbour edges for edges" + Mesh3_edge_links[...] = [ + [1, 2, 3, -99, -99], + [0, 2, 5, 8, -99], + [3, 0, 1, 5, 8], + [4, 2, 0, -99, -99], + [ + 3, + 5, + 6, + -99, + -99, + ], + [4, 6, 2, 1, 8], + [ + 4, + 5, + 7, + -99, + -99, + ], + [ + 6, + 8, + -99, + -99, + -99, + ], + [7, 5, 2, 1, -99], + ] + + n.close() + return filename + + def _make_aggregation_value(filename): """Create an aggregation variable with 'unique_values'.""" n = netCDF4.Dataset(filename, "w") @@ -2342,6 +2486,7 @@ def _make_aggregation_value(filename): ugrid_1 = _make_ugrid_1("ugrid_1.nc") ugrid_2 = _make_ugrid_2("ugrid_2.nc") +ugrid_3 = _make_ugrid_3("ugrid_3.nc") aggregation_value = _make_aggregation_value("aggregation_value.nc") diff --git a/cfdm/test/test_DomainTopology.py b/cfdm/test/test_DomainTopology.py index 27135c276..8a6d3632f 100644 --- a/cfdm/test/test_DomainTopology.py +++ b/cfdm/test/test_DomainTopology.py @@ -8,23 +8,13 @@ import cfdm -# Create test domain topology object -c = cfdm.DomainTopology() -c.set_properties({"long_name": "Maps every face to its corner nodes"}) -c.nc_set_variable("Mesh2_face_nodes") -data = cfdm.Data( - [[2, 3, 1, 0], [6, 7, 3, 2], [1, 3, 8, -99]], - dtype="i4", -) -data.masked_values(-99, inplace=True) -c.set_data(data) -c.set_cell("face") - class DomainTopologyTest(unittest.TestCase): """Unit test for the DomainTopology class.""" - d = c + face = cfdm.example_field(8).domain_topology() + edge = cfdm.example_field(9).domain_topology() + point = cfdm.example_field(10).domain_topology() def setUp(self): """Preparations called immediately before each test method.""" @@ -41,7 +31,7 @@ def setUp(self): def test_DomainTopology__repr__str__dump(self): """Test all means of DomainTopology inspection.""" - d = self.d + d = self.face self.assertEqual(repr(d), "") self.assertEqual(str(d), "cell:face(3, 4) ") self.assertEqual( @@ -53,17 +43,17 @@ def test_DomainTopology__repr__str__dump(self): def test_DomainTopology_copy(self): """Test the copy of DomainTopology.""" - d = self.d + d = self.face self.assertTrue(d.equals(d.copy())) def test_DomainTopology_data(self): """Test the data of DomainTopology.""" - d = self.d + d = self.face self.assertEqual(d.ndim, 1) def test_DomainTopology_cell(self): """Test the 'cell' methods of DomainTopology.""" - d = self.d.copy() + d = self.face.copy() self.assertTrue(d.has_cell()) self.assertEqual(d.get_cell(), "face") self.assertEqual(d.del_cell(), "face") @@ -83,7 +73,7 @@ def test_DomainTopology_cell(self): def test_DomainTopology_transpose(self): """Test the 'transpose' method of DomainTopology.""" - d = self.d.copy() + d = self.face.copy() e = d.transpose() self.assertTrue(d.equals(e)) self.assertIsNone(d.transpose(inplace=True)) @@ -199,6 +189,85 @@ def test_DomainTopology_normalise(self): self.assertTrue((d0.mask == n.mask).all()) self.assertTrue((d0 == n).all()) + def test_DomainTopology_sort(self): + """Test cfdm.DomainTopology.sort.""" + # Test sorting of "point" + d = self.point + e = d.sort() + self.assertIsInstance(e, cfdm.DomainTopology) + + result = cfdm.Data( + [ + [0, 1, 2, -99, -99], + [1, 0, 3, 6, -99], + [2, 0, 3, 4, -99], + [3, 1, 2, 5, 6], + [4, 2, 5, -99, -99], + [5, 3, 4, -99, -99], + [6, 1, 3, -99, -99], + ], + dtype="i4", + mask_value=-99, + ) + self.assertTrue(e.data.equals(result)) + + e = d[::-1] + self.assertFalse(e.equals(d)) + e = e.sort() + self.assertTrue(e.data.equals(result)) + + # Test sorting of "edge" + d = self.edge + e = d.sort() + self.assertIsInstance(e, cfdm.DomainTopology) + + result = cfdm.Data( + [ + [0, 1], + [0, 2], + [1, 3], + [1, 6], + [2, 3], + [2, 4], + [3, 5], + [3, 6], + [4, 5], + ], + dtype="i4", + mask_value=-99, + ) + self.assertTrue(e.data.equals(result)) + + e = d[::-1] + self.assertFalse(e.equals(d)) + e = e.sort() + self.assertTrue(e.data.equals(result)) + + # Can't sort "face" + with self.assertRaises(ValueError): + self.face.sort() + + def test_DomainTopology_to_edge(self): + """Test cfdm.DomainTopology.to_edge.""" + # Test "edge" to_edge + d = self.edge + e = d.to_edge() + self.assertIsInstance(e, cfdm.DomainTopology) + self.assertTrue(e.equals(d)) + self.assertTrue(d.to_edge(sort=True).equals(d.sort())) + + # Test "point" to_edge + d = self.point + e = d.to_edge(sort=True) + self.assertIsInstance(e, cfdm.DomainTopology) + self.assertTrue(e.data.equals(self.edge.sort().data)) + + # Test "face" to_edge + d = self.face + e = d.to_edge(sort=True) + self.assertIsInstance(e, cfdm.DomainTopology) + self.assertTrue(e.data.equals(self.edge.sort().data)) + if __name__ == "__main__": print("Run date:", datetime.datetime.now()) diff --git a/cfdm/test/test_UGRID.py b/cfdm/test/test_UGRID.py index eac354459..7c8b6112f 100644 --- a/cfdm/test/test_UGRID.py +++ b/cfdm/test/test_UGRID.py @@ -1,10 +1,12 @@ import atexit import datetime import faulthandler +import itertools import os import tempfile import unittest +import netCDF4 import numpy as np faulthandler.enable() # to debug seg faults and timeouts @@ -14,12 +16,12 @@ warnings = False # Set up temporary files -n_tmpfiles = 1 +n_tmpfiles = 2 tmpfiles = [ - tempfile.mkstemp("_test_read_write.nc", dir=os.getcwd())[1] + tempfile.mkstemp("_test_ugrid.nc", dir=os.getcwd())[1] for i in range(n_tmpfiles) ] -[tmpfile1] = tmpfiles +[tmpfile, tmpfile1] = tmpfiles def _remove_tmpfiles(): @@ -34,6 +36,31 @@ def _remove_tmpfiles(): atexit.register(_remove_tmpfiles) +def n_mesh_variables(filename): + """Return the number of mesh variables in the file.""" + nc = netCDF4.Dataset(filename, "r") + n = 0 + for v in nc.variables.values(): + try: + v.getncattr("topology_dimension") + except AttributeError: + pass + else: + n += 1 + + nc.close() + return n + + +def combinations(face, edge, point): + """Return the combination for field/domain indexing.""" + return [ + i + for n in range(1, 4) + for i in itertools.permutations([face, edge, point], n) + ] + + class UGRIDTest(unittest.TestCase): """Test UGRID field constructs.""" @@ -45,6 +72,10 @@ class UGRIDTest(unittest.TestCase): os.path.dirname(os.path.abspath(__file__)), "ugrid_2.nc" ) + filename3 = os.path.join( + os.path.dirname(os.path.abspath(__file__)), "ugrid_3.nc" + ) + def setUp(self): """Preparations called immediately before each test method.""" # Disable log messages to silence expected warnings @@ -76,10 +107,6 @@ def test_UGRID_read(self): g.cell_connectivity().get_connectivity(), "edge" ) - # Check that all fields have the same mesh id - mesh_ids1 = set(g.get_mesh_id() for g in f1) - self.assertEqual(len(mesh_ids1), 1) - f2 = cfdm.read(self.filename2) self.assertEqual(len(f2), 3) for g in f2: @@ -98,13 +125,6 @@ def test_UGRID_read(self): g.cell_connectivity().get_connectivity(), "edge" ) - # Check that all fields have the same mesh id - mesh_ids2 = set(g.get_mesh_id() for g in f2) - self.assertEqual(len(mesh_ids2), 1) - - # Check that the different files have different mesh ids - self.assertNotEqual(mesh_ids1, mesh_ids2) - def test_UGRID_data(self): """Test reading of UGRID data.""" node1, face1, edge1 = cfdm.read(self.filename1) @@ -177,9 +197,102 @@ def test_read_UGRID_domain(self): g.cell_connectivity().get_connectivity(), "edge" ) - # Check that all domains have the same mesh id - mesh_ids1 = set(g.get_mesh_id() for g in d1) - self.assertEqual(len(mesh_ids1), 1) + def test_read_write_UGRID_field(self): + """Test the cfdm.read and cfdm.write with UGRID fields.""" + # Face, edge, and point fields that are all part of the same + # UGRID mesh + ugrid = cfdm.example_fields(8, 9, 10) + face, edge, point = (0, 1, 2) + + # Test for equality with the fields defined in memory. Only + # works for face and edge fields. + for cell in (face, edge): + f = ugrid[cell] + cfdm.write(f, tmpfile) + g = cfdm.read(tmpfile) + self.assertEqual(len(g), 1) + self.assertTrue(g[0].equals(f)) + + # Test round-tripping of field combinations + for cells in combinations(face, edge, point): + f = [] + for cell in cells: + f.append(ugrid[cell]) + + cfdm.write(f, tmpfile) + + # Check that there's only one mesh variable in the file + self.assertEqual(n_mesh_variables(tmpfile), 1) + + g = cfdm.read(tmpfile) + self.assertEqual(len(g), len(f)) + + cfdm.write(g, tmpfile1) + + # Check that there's only one mesh variable in the file + self.assertEqual(n_mesh_variables(tmpfile1), 1) + + h = cfdm.read(tmpfile1) + self.assertEqual(len(h), len(g)) + self.assertTrue(h[0].equals(g[0])) + + def test_read_write_UGRID_domain(self): + """Test the cfdm.read and cfdm.write with UGRID domains.""" + # Face, edge, and point fields/domains that are all part of + # the same UGRID mesh + ugrid = [f.domain for f in cfdm.example_fields(8, 9, 10)] + face, edge, point = (0, 1, 2) + + # Test for equality with the fields defined in memory. Only + # works for face and edge domains. + for cell in (face, edge): + d = ugrid[cell] + cfdm.write(d, tmpfile) + e = cfdm.read(tmpfile, domain=True) + self.assertEqual(len(e), 2) + self.assertTrue(e[0].equals(d)) + self.assertEqual(e[1].domain_topology().get_cell(), "point") + + # Test round-tripping of domain combinations for the + # example_field domains, and also the domain read from + # 'ugrid_3.nc'. + for iteration in ("memory", "file"): + for cells in combinations(face, edge, point): + d = [] + for cell in cells: + d.append(ugrid[cell]) + + if point not in cells: + # When we write a non-point domains, we also get + # the point locations. + d.append(ugrid[point]) + elif cells == (point,): + # When we write a point domain on its own, we also + # get the edge location. + d.append(ugrid[edge]) + + cfdm.write(d, tmpfile) + + # Check that there's only one mesh variable in the file + self.assertEqual(n_mesh_variables(tmpfile), 1) + + e = cfdm.read(tmpfile, domain=True) + + self.assertEqual(len(e), len(d)) + + cfdm.write(e, tmpfile1) + + # Check that there's only one mesh variable in the file + self.assertEqual(n_mesh_variables(tmpfile1), 1) + + f = cfdm.read(tmpfile1, domain=True) + self.assertEqual(len(f), len(e)) + for i, j in zip(f, e): + self.assertTrue(i.equals(j)) + + # Set up for the 'file' iteration + ugrid = cfdm.read(self.filename3, domain=True) + face, edge, point = (2, 1, 0) if __name__ == "__main__": diff --git a/cfdm/test/test_read_write.py b/cfdm/test/test_read_write.py index 14c55d6d1..67cfda508 100644 --- a/cfdm/test/test_read_write.py +++ b/cfdm/test/test_read_write.py @@ -273,9 +273,8 @@ def test_write_netcdf_mode(self): if fmt == "NETCDF4_CLASSIC" and ex_field_n in (6, 7): continue - print( - "TODOUGRID: excluding example fields 8, 9, 10 until writing UGRID is enabled" - ) + # Exclude UGRID fields, as we deal with them in + # test_UGRID.py if ex_field_n in (8, 9, 10): continue @@ -310,9 +309,9 @@ def test_write_netcdf_mode(self): # Now do the same test, but appending all of the example fields in # one operation rather than one at a time, to check that it works. cfdm.write(g, tmpfile, fmt=fmt, mode="w") # 1. overwrite to wipe - print( - "TODOUGRID: excluding example fields 8, 9, 10 until writing UGRID is enabled" - ) + + # Exclude UGRID fields, as we deal with them in + # test_UGRID.py append_ex_fields = cfdm.example_fields(*range(8)) del append_ex_fields[1] # note: can remove after Issue #141 closed # Note: can remove this del when Issue #140 is closed: diff --git a/cfdm/test/test_zarr.py b/cfdm/test/test_zarr.py index 3f24ea303..5943ff9a6 100644 --- a/cfdm/test/test_zarr.py +++ b/cfdm/test/test_zarr.py @@ -161,10 +161,10 @@ def test_zarr_write_append(self): """Test in append mode with Zarr.""" # Check that append mode does not work for Zarr f = self.f0 - cfdm.write(f, tmpdir1, fmt='ZARR3') + cfdm.write(f, tmpdir1, fmt="ZARR3") with self.assertRaises(ValueError): - cfdm.write(f, tmpdir1, fmt='ZARR3', mode="a") - + cfdm.write(f, tmpdir1, fmt="ZARR3", mode="a") + def test_zarr_groups_1(self): """Test for the general handling of Zarr hierarchical groups.""" f = cfdm.example_field(1)