Ugrid save

lib/iris/_representation/cube_summary.py

@@ -259,8 +259,17 @@ def __init__(self, cube, shorten=False, name_padding=35):
         vector_dim_coords = [
             coord for coord in dim_coords if id(coord) not in scalar_coord_ids
         ]
+        if cube.mesh is None:
+            mesh_coords = []
+        else:
+            mesh_coords = [
+                coord for coord in aux_coords if hasattr(coord, "mesh")
+            ]
+
         vector_aux_coords = [
-            coord for coord in aux_coords if id(coord) not in scalar_coord_ids
+            coord
+            for coord in aux_coords
+            if (id(coord) not in scalar_coord_ids and coord not in mesh_coords)
         ]
         vector_derived_coords = [
             coord
@@ -300,6 +309,7 @@ def add_vector_section(title, contents, iscoord=True):
             )
 
         add_vector_section("Dimension coordinates:", vector_dim_coords)
+        add_vector_section("Mesh coordinates:", mesh_coords)
         add_vector_section("Auxiliary coordinates:", vector_aux_coords)
         add_vector_section("Derived coordinates:", vector_derived_coords)
         add_vector_section("Cell measures:", vector_cell_measures, False)

lib/iris/experimental/ugrid/__init__.py

@@ -522,6 +522,8 @@ def __eq__(self, other):
             if hasattr(other, "metadata"):
                 # metadata comparison
                 eq = self.metadata == other.metadata
+                if eq:
+                    eq = self.shape == other.shape
                 if eq:
                     eq = (
                         self.indices_by_src() == other.indices_by_src()

lib/iris/tests/integration/experimental/test_ugrid_save.py

@@ -0,0 +1,127 @@
+# Copyright Iris contributors
+#
+# This file is part of Iris and is released under the LGPL license.
+# See COPYING and COPYING.LESSER in the root of the repository for full
+# licensing details.
+"""
+Integration tests for NetCDF-UGRID file saving.
+
+"""
+# Import iris.tests first so that some things can be initialised before
+# importing anything else.
+import iris.tests as tests  # isort:skip
+
+import glob
+from pathlib import Path
+import shutil
+from subprocess import check_call
+import tempfile
+
+import iris
+from iris.experimental.ugrid import PARSE_UGRID_ON_LOAD
+import iris.fileformats.netcdf
+from iris.tests import IrisTest
+from iris.tests.stock.netcdf import _add_standard_data
+
+
+class TestBasicSave(IrisTest):
+    @classmethod
+    def setUpClass(cls):
+        cls.temp_dir = Path(tempfile.mkdtemp())
+        cls.examples_dir = (
+            Path(__file__).absolute().parent / "ugrid_conventions_examples"
+        )
+        example_paths = glob.glob(str(cls.examples_dir / "*ex*.cdl"))
+        example_names = [
+            str(Path(filepath).name).split("_")[1]  # = "ex<N>"
+            for filepath in example_paths
+        ]
+        cls.example_names_paths = {
+            name: path for name, path in zip(example_names, example_paths)
+        }
+
+    @classmethod
+    def tearDownClass(cls):
+        shutil.rmtree(cls.temp_dir)
+
+    def test_example_result_cdls(self):
+        # Snapshot the result of saving the example cases.
+        for ex_name, filepath in self.example_names_paths.items():
+            target_ncfile_path = str(self.temp_dir / f"{ex_name}.nc")
+            # Create a netcdf file from the test CDL.
+            check_call(
+                f"ncgen {filepath} -k4 -o {target_ncfile_path}", shell=True
+            )
+            # Fill in blank data-variables.
+            _add_standard_data(target_ncfile_path)
+            # Load as Iris data
+            with PARSE_UGRID_ON_LOAD.context():
+                cubes = iris.load(target_ncfile_path)
+            # Re-save, to check the save behaviour.
+            resave_ncfile_path = str(self.temp_dir / f"{ex_name}_resaved.nc")
+            iris.save(cubes, resave_ncfile_path)
+            # Check the output against a CDL snapshot.
+            refdir_relpath = (
+                "integration/experimental/ugrid_save/TestBasicSave/"
+            )
+            reffile_name = str(Path(filepath).name).replace(".nc", ".cdl")
+            reffile_path = refdir_relpath + reffile_name
+            self.assertCDL(resave_ncfile_path, reference_filename=reffile_path)
+
+    def test_example_roundtrips(self):
+        # Check that save-and-loadback leaves Iris data unchanged,
+        # for data derived from each UGRID example CDL.
+        for ex_name, filepath in self.example_names_paths.items():
+            target_ncfile_path = str(self.temp_dir / f"{ex_name}.nc")
+            # Create a netcdf file from the test CDL.
+            check_call(
+                f"ncgen {filepath} -k4 -o {target_ncfile_path}", shell=True
+            )
+            # Fill in blank data-variables.
+            _add_standard_data(target_ncfile_path)
+            # Load the original as Iris data
+            with PARSE_UGRID_ON_LOAD.context():
+                orig_cubes = iris.load(target_ncfile_path)
+
+            if "ex4" in ex_name:
+                # Discard the extra formula terms component cubes
+                # Saving these does not do what you expect
+                orig_cubes = orig_cubes.extract("datavar")
+
+            # Save-and-load-back to compare the Iris saved result.
+            resave_ncfile_path = str(self.temp_dir / f"{ex_name}_resaved.nc")
+            iris.save(orig_cubes, resave_ncfile_path)
+            with PARSE_UGRID_ON_LOAD.context():
+                savedloaded_cubes = iris.load(resave_ncfile_path)
+
+            # This should match the original exactly
+            # ..EXCEPT for our inability to compare meshes.
+            for orig, reloaded in zip(orig_cubes, savedloaded_cubes):
+                for cube in (orig, reloaded):
+                    # Remove conventions attributes, which may differ.
+                    cube.attributes.pop("Conventions", None)
+                    # Remove var-names, which may differ.
+                    cube.var_name = None
+
+                # Compare the mesh contents (as we can't compare actual meshes)
+                self.assertEqual(orig.location, reloaded.location)
+                orig_mesh = orig.mesh
+                reloaded_mesh = reloaded.mesh
+                self.assertEqual(
+                    orig_mesh.all_coords, reloaded_mesh.all_coords
+                )
+                self.assertEqual(
+                    orig_mesh.all_connectivities,
+                    reloaded_mesh.all_connectivities,
+                )
+                # Index the cubes to replace meshes with meshcoord-derived aux coords.
+                # This needs [:0] on the mesh dim, so do that on all dims.
+                keys = tuple([slice(0, None)] * orig.ndim)
+                orig = orig[keys]
+                reloaded = reloaded[keys]
+                # Resulting cubes, with collapsed mesh, should be IDENTICAL.
+                self.assertEqual(orig, reloaded)
+
+
+if __name__ == "__main__":
+    tests.main()

lib/iris/tests/integration/experimental/ugrid_conventions_examples/README.txt

@@ -0,0 +1,16 @@
+Examples generated from CDL example sections in UGRID conventions v1.0 
+  ( see webpage: https://ugrid-conventions.github.io/ugrid-conventions/ )
+
+CHANGES:
+    * added a data-var to all examples, for ease of iris-roundtripping
+    * EX4 :
+        - had a couple of missing ";"s at lineends
+        - the formula terms (depth+surface) should map to 'Mesh2_layers', and not to the mesh at all.
+            - use Mesh2d_layers dim, and have no 'mesh' or 'location'
+    * "EX4a" -- possibly (future) closer mix of hybrid-vertical and mesh dimensions
+        - *don't* think we can have a hybrid coord ON the mesh dimension
+            - mesh being a vertical location (only) seems to make no sense
+            - .. and implies that the mesh is 1d and ordered, which is not really unstructured at all
+        - *could* have hybrid-height with the _orography_ mapping to the mesh
+            - doesn't match the UGRID examples, but see : iris.tests.unit.fileformats.netcdf.test_Saver__ugrid.TestSaveUgrid__cube.test_nonmesh_hybrid_dim
+

lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex1_1d_mesh.cdl

@@ -0,0 +1,55 @@
+netcdf ex1_1d_mesh {
+dimensions:
+nMesh1_node = 5 ; // nNodes
+nMesh1_edge = 4 ; // nEdges
+
+Two = 2;
+
+variables:
+// Mesh topology
+integer Mesh1 ;
+Mesh1:cf_role = "mesh_topology" ;
+Mesh1:long_name = "Topology data of 1D network" ;
+Mesh1:topology_dimension = 1 ;
+Mesh1:node_coordinates = "Mesh1_node_x Mesh1_node_y" ;
+Mesh1:edge_node_connectivity = "Mesh1_edge_nodes" ;
+Mesh1:edge_coordinates = "Mesh1_edge_x Mesh1_edge_y" ; // optional attribute
+integer Mesh1_edge_nodes(nMesh1_edge, Two) ;
+Mesh1_edge_nodes:cf_role = "edge_node_connectivity" ;
+Mesh1_edge_nodes:long_name = "Maps every edge/link to the two nodes that it connects." ;
+Mesh1_edge_nodes:start_index = 1 ;
+
+// Mesh node coordinates
+double Mesh1_node_x(nMesh1_node) ;
+Mesh1_node_x:standard_name = "longitude" ;
+Mesh1_node_x:long_name = "Longitude of 1D network nodes." ;
+Mesh1_node_x:units = "degrees_east" ;
+double Mesh1_node_y(nMesh1_node) ;
+Mesh1_node_y:standard_name = "latitude" ;
+Mesh1_node_y:long_name = "Latitude of 1D network nodes." ;
+Mesh1_node_y:units = "degrees_north" ;
+
+// Optional mesh edge coordinate variables
+double Mesh1_edge_x(nMesh1_edge) ;
+Mesh1_edge_x:standard_name = "longitude" ;
+Mesh1_edge_x:long_name = "Characteristic longitude of 1D network edge (e.g. midpoint of the edge)." ;
+Mesh1_edge_x:units = "degrees_east" ;
+Mesh1_edge_x:bounds = "Mesh1_edge_xbnds" ;
+double Mesh1_edge_y(nMesh1_edge) ;
+Mesh1_edge_y:standard_name = "latitude" ;
+Mesh1_edge_y:long_name = "Characteristic latitude of 1D network edge (e.g. midpoint of the edge)." ;
+Mesh1_edge_y:units = "degrees_north" ;
+Mesh1_edge_y:bounds = "Mesh1_edge_ybnds" ;
+double Mesh1_edge_xbnds(nMesh1_edge,Two) ;
+Mesh1_edge_xbnds:standard_name = "longitude" ;
+Mesh1_edge_xbnds:long_name = "Longitude bounds of 1D network edge (i.e. begin and end longitude)." ;
+Mesh1_edge_xbnds:units = "degrees_east" ;
+double Mesh1_edge_ybnds(nMesh1_edge,Two) ;
+Mesh1_edge_ybnds:standard_name = "latitude" ;
+Mesh1_edge_ybnds:long_name = "Latitude bounds of 1D network edge (i.e. begin and end latitude)." ;
+Mesh1_edge_ybnds:units = "degrees_north" ;
+
+float datavar(nMesh1_edge) ;
+  datavar:mesh = "Mesh1" ;
+  datavar:location = "edge" ;
+}

lib/iris/tests/integration/experimental/ugrid_conventions_examples/ugrid_ex2_2d_triangular.cdl

@@ -0,0 +1,84 @@
+netcdf ex2_2d_triangular {
+dimensions:
+nMesh2_node = 4 ; // nNodes
+nMesh2_edge = 5 ; // nEdges
+nMesh2_face = 2 ; // nFaces
+
+Two = 2 ;
+Three = 3 ;
+
+variables:
+// Mesh topology
+integer Mesh2 ;
+Mesh2:cf_role = "mesh_topology" ;
+Mesh2:long_name = "Topology data of 2D unstructured mesh" ;
+Mesh2:topology_dimension = 2 ;
+Mesh2:node_coordinates = "Mesh2_node_x Mesh2_node_y" ;
+Mesh2:face_node_connectivity = "Mesh2_face_nodes" ;
+Mesh2:face_dimension = "nMesh2_face" ;
+Mesh2:edge_node_connectivity = "Mesh2_edge_nodes" ; // attribute required if variables will be defined on edges
+Mesh2:edge_dimension = "nMesh2_edge" ;
+Mesh2:edge_coordinates = "Mesh2_edge_x Mesh2_edge_y" ; // optional attribute (requires edge_node_connectivity)
+Mesh2:face_coordinates = "Mesh2_face_x Mesh2_face_y" ; // optional attribute
+Mesh2:face_edge_connectivity = "Mesh2_face_edges" ; // optional attribute (requires edge_node_connectivity)
+Mesh2:face_face_connectivity = "Mesh2_face_links" ; // optional attribute
+Mesh2:edge_face_connectivity = "Mesh2_edge_face_links" ; // optional attribute (requires edge_node_connectivity)
+integer Mesh2_face_nodes(nMesh2_face, Three) ;
+Mesh2_face_nodes:cf_role = "face_node_connectivity" ;
+Mesh2_face_nodes:long_name = "Maps every triangular face to its three corner nodes." ;
+Mesh2_face_nodes:start_index = 1 ;
+integer Mesh2_edge_nodes(nMesh2_edge, Two) ;
+Mesh2_edge_nodes:cf_role = "edge_node_connectivity" ;
+Mesh2_edge_nodes:long_name = "Maps every edge to the two nodes that it connects." ;
+Mesh2_edge_nodes:start_index = 1 ;
+
+// Optional mesh topology variables
+integer Mesh2_face_edges(nMesh2_face, Three) ;
+Mesh2_face_edges:cf_role = "face_edge_connectivity" ;
+Mesh2_face_edges:long_name = "Maps every triangular face to its three edges." ;
+Mesh2_face_edges:start_index = 1 ;
+integer Mesh2_face_links(nMesh2_face, Three) ;
+Mesh2_face_links:cf_role = "face_face_connectivity" ;
+Mesh2_face_links:long_name = "neighbor faces for faces" ;
+Mesh2_face_links:start_index = 1 ;
+Mesh2_face_links:_FillValue = -999 ;
+Mesh2_face_links:comment = "missing neighbor faces are indicated using _FillValue" ;
+integer Mesh2_edge_face_links(nMesh2_edge, Two) ;
+Mesh2_edge_face_links:cf_role = "edge_face_connectivity" ;
+Mesh2_edge_face_links:long_name = "neighbor faces for edges" ;
+Mesh2_edge_face_links:start_index = 1 ;
+Mesh2_edge_face_links:_FillValue = -999 ;
+Mesh2_edge_face_links:comment = "missing neighbor faces are indicated using _FillValue" ;
+
+// Mesh node coordinates
+double Mesh2_node_x(nMesh2_node) ;
+Mesh2_node_x:standard_name = "longitude" ;
+Mesh2_node_x:long_name = "Longitude of 2D mesh nodes." ;
+Mesh2_node_x:units = "degrees_east" ;
+double Mesh2_node_y(nMesh2_node) ;
+Mesh2_node_y:standard_name = "latitude" ;
+Mesh2_node_y:long_name = "Latitude of 2D mesh nodes." ;
+Mesh2_node_y:units = "degrees_north" ;
+
+// Optional mesh face and edge coordinate variables
+double Mesh2_face_x(nMesh2_face) ;
+Mesh2_face_x:standard_name = "longitude" ;
+Mesh2_face_x:long_name = "Characteristics longitude of 2D mesh triangle (e.g. circumcenter coordinate)." ;
+Mesh2_face_x:units = "degrees_east" ;
+double Mesh2_face_y(nMesh2_face) ;
+Mesh2_face_y:standard_name = "latitude" ;
+Mesh2_face_y:long_name = "Characteristics latitude of 2D mesh triangle (e.g. circumcenter coordinate)." ;
+Mesh2_face_y:units = "degrees_north" ;
+double Mesh2_edge_x(nMesh2_edge) ;
+Mesh2_edge_x:standard_name = "longitude" ;
+Mesh2_edge_x:long_name = "Characteristic longitude of 2D mesh edge (e.g. midpoint of the edge)." ;
+Mesh2_edge_x:units = "degrees_east" ;
+double Mesh2_edge_y(nMesh2_edge) ;
+Mesh2_edge_y:standard_name = "latitude" ;
+Mesh2_edge_y:long_name = "Characteristic latitude of 2D mesh edge (e.g. midpoint of the edge)." ;
+Mesh2_edge_y:units = "degrees_north" ;
+
+float datavar(nMesh2_face) ;
+  datavar:mesh = "Mesh2" ;
+  datavar:location = "face" ;
+}