SciTools · ESadek-MO · Jun 20, 2024 · Jun 12, 2024 · Jun 12, 2024 · Jun 12, 2024
diff --git a/docs/src/whatsnew/latest.rst b/docs/src/whatsnew/latest.rst
@@ -44,6 +44,10 @@ This document explains the changes made to Iris for this release
    :func:`~iris.cube.Cube.rolling_window`.  (:issue:`5777`, :pull:`5825`)
 
 
+#. `@pp-mo`_ corrected the use of mesh dimensions when saving with multiple
+   meshes.  (:issue:`5908`, :pull:`6004`)
+
+
 💣 Incompatible Changes
 =======================
 

diff --git a/lib/iris/experimental/ugrid/load.py b/lib/iris/experimental/ugrid/load.py
@@ -485,9 +485,27 @@ def _build_mesh_coords(mesh, cf_var):
         "edge": mesh.edge_dimension,
         "face": mesh.face_dimension,
     }
-    mesh_dim_name = element_dimensions[cf_var.location]
-    # (Only expecting 1 mesh dimension per cf_var).
-    mesh_dim = cf_var.dimensions.index(mesh_dim_name)
+    location = getattr(cf_var, "location", "<empty>")
+    if location is None or location not in element_dimensions:
+        # We should probably issue warnings and recover, but that is too much
+        # work.  Raising a more intelligible error is easy to do though.
+        msg = (
+            f"mesh data variable {cf_var.name!r} has an invalid "
+            f"location={location!r}."
+        )
+        raise ValueError(msg)
+    mesh_dim_name = element_dimensions.get(location)
+    if mesh_dim_name is None:
+        msg = f"mesh {mesh.name!r} has no {location} dimension."
+        raise ValueError(msg)
+    if mesh_dim_name in cf_var.dimensions:
+        mesh_dim = cf_var.dimensions.index(mesh_dim_name)
+    else:
+        msg = (
+            f"mesh data variable {cf_var.name!r} does not have the "
+            f"{location} mesh dimension {mesh_dim_name!r}, in its dimensions."
+        )
+        raise ValueError(msg)
 
     mesh_coords = mesh.to_MeshCoords(location=cf_var.location)
     return mesh_coords, mesh_dim
diff --git a/lib/iris/fileformats/netcdf/saver.py b/lib/iris/fileformats/netcdf/saver.py
@@ -1216,13 +1216,6 @@ def record_dimension(names_list, dim_name, length, matching_coords=None):
                     assert len(dim_coords) == 1
                     dim_element = dim_coords[0]
                     dim_name = self._dim_names_and_coords.name(dim_element)
-                    if dim_name is not None:
-                        # For mesh-identifying coords, we require the *same*
-                        # coord, not an identical one (i.e. "is" not "==")
-                        stored_coord = self._dim_names_and_coords.coord(dim_name)
-                        if dim_element is not stored_coord:
-                            # This is *not* a proper match after all.
-                            dim_name = None
                     if dim_name is None:
                         # No existing dim matches this, so assign a new name
                         if location == "node":

diff --git a/lib/iris/tests/unit/experimental/ugrid/load/test_load_meshes.py b/lib/iris/tests/unit/experimental/ugrid/load/test_load_meshes.py
@@ -27,50 +27,58 @@ def tearDownModule():
         rmtree(TMP_DIR)
 
 
-def cdl_to_nc(cdl):
-    cdl_path = str(TMP_DIR / "tst.cdl")
-    nc_path = str(TMP_DIR / f"{uuid4()}.nc")
+def cdl_to_nc(cdl, tmpdir=None):
+    if tmpdir is None:
+        tmpdir = TMP_DIR
+    cdl_path = str(tmpdir / "tst.cdl")
+    nc_path = str(tmpdir / f"{uuid4()}.nc")
     # Use ncgen to convert this into an actual (temporary) netCDF file.
     ncgen_from_cdl(cdl_str=cdl, cdl_path=cdl_path, nc_path=nc_path)
     return nc_path
 
 
-class TestsBasic(tests.IrisTest):
+_TEST_CDL_HEAD = """
+netcdf mesh_test {
+    dimensions:
+        node = 3 ;
+        face = 1 ;
+        vertex = 3 ;
+        levels = 2 ;
+    variables:
+        int mesh ;
+            mesh:cf_role = "mesh_topology" ;
+            mesh:topology_dimension = 2 ;
+            mesh:node_coordinates = "node_x node_y" ;
+            mesh:face_node_connectivity = "face_nodes" ;
+        float node_x(node) ;
+            node_x:standard_name = "longitude" ;
+        float node_y(node) ;
+            node_y:standard_name = "latitude" ;
+        int face_nodes(face, vertex) ;
+            face_nodes:cf_role = "face_node_connectivity" ;
+            face_nodes:start_index = 0 ;
+        int levels(levels) ;
+        float node_data(levels, node) ;
+            node_data:coordinates = "node_x node_y" ;
+            node_data:location = "node" ;
+            node_data:mesh = "mesh" ;
+"""
+
+_TEST_CDL_TAIL = """
+data:
+        mesh = 0;
+        node_x = 0., 2., 1.;
+        node_y = 0., 0., 1.;
+        face_nodes = 0, 1, 2;
+        levels = 1, 2;
+        node_data = 0., 0., 0.;
+    }
+"""
+
+
+class TestLoadErrors(tests.IrisTest):
     def setUp(self):
-        self.ref_cdl = """
-            netcdf mesh_test {
-                dimensions:
-                    node = 3 ;
-                    face = 1 ;
-                    vertex = 3 ;
-                    levels = 2 ;
-                variables:
-                    int mesh ;
-                        mesh:cf_role = "mesh_topology" ;
-                        mesh:topology_dimension = 2 ;
-                        mesh:node_coordinates = "node_x node_y" ;
-                        mesh:face_node_connectivity = "face_nodes" ;
-                    float node_x(node) ;
-                        node_x:standard_name = "longitude" ;
-                    float node_y(node) ;
-                        node_y:standard_name = "latitude" ;
-                    int face_nodes(face, vertex) ;
-                        face_nodes:cf_role = "face_node_connectivity" ;
-                        face_nodes:start_index = 0 ;
-                    int levels(levels) ;
-                    float node_data(levels, node) ;
-                        node_data:coordinates = "node_x node_y" ;
-                        node_data:location = "node" ;
-                        node_data:mesh = "mesh" ;
-                data:
-                    mesh = 0;
-                    node_x = 0., 2., 1.;
-                    node_y = 0., 0., 1.;
-                    face_nodes = 0, 1, 2;
-                    levels = 1, 2;
-                    node_data = 0., 0., 0.;
-                }
-            """
+        self.ref_cdl = _TEST_CDL_HEAD + _TEST_CDL_TAIL
         self.nc_path = cdl_to_nc(self.ref_cdl)
 
     def add_second_mesh(self):

diff --git a/lib/iris/tests/unit/experimental/ugrid/load/test_meshload_checks.py b/lib/iris/tests/unit/experimental/ugrid/load/test_meshload_checks.py
@@ -0,0 +1,83 @@
+# Copyright Iris contributors
+#
+# This file is part of Iris and is released under the BSD license.
+# See LICENSE in the root of the repository for full licensing details.
+"""Unit tests for mesh handling within iris netcdf loads."""
+
+import pytest
+
+import iris
+from iris.experimental.ugrid.load import PARSE_UGRID_ON_LOAD
+
+from .test_load_meshes import (
+    _TEST_CDL_HEAD,
+    _TEST_CDL_TAIL,
+    cdl_to_nc,
+)
+
+
+class TestMeshLoad:
+    def _create_testnc(self, location="node", meshdim="node"):
+        # Add an extra (possibly mal-formed) mesh data to the testfile.
+        if location is None:
+            location_cdl = ""
+        else:
+            location_cdl = f'extra_data:location = "{location}" ;'
+
+        extra_cdl = f"""
+        float extra_data(levels, {meshdim}) ;
+            extra_data:coordinates = "node_x node_y" ;
+            {location_cdl}
+            extra_data:mesh = "mesh" ;
+        """
+        # Insert this into the definitions part of the 'standard' testfile CDL
+        extended_cdl = _TEST_CDL_HEAD + extra_cdl + _TEST_CDL_TAIL
+        testfile_path = cdl_to_nc(extended_cdl, tmpdir=self.tmpdir)
+        return testfile_path
+
+    @pytest.fixture(params=["nolocation", "badlocation", "baddim"])
+    def failnc(self, request, tmp_path_factory):
+        self.param = request.param
+        kwargs = {}
+        if self.param == "nolocation":
+            kwargs["location"] = None
+        elif self.param == "badlocation":
+            kwargs["location"] = "invalid_location"
+        elif self.param == "baddim":
+            kwargs["meshdim"] = "vertex"
+        else:
+            raise ValueError(f"unexpected param: {self.param}")
+
+        self.tmpdir = tmp_path_factory.mktemp("meshload")
+        return self._create_testnc(**kwargs)
+
+    def test_extrameshvar__ok(self, tmp_path_factory):
+        # Check that the default cdl construction loads OK
+        self.tmpdir = tmp_path_factory.mktemp("meshload")
+        testnc = self._create_testnc()
+        with PARSE_UGRID_ON_LOAD.context():
+            iris.load(testnc)
+
+    def test_extrameshvar__fail(self, failnc):
+        # Check that the expected errors are raised in various cases.
+        param = self.param
+        if param == "nolocation":
+            match_msg = (
+                "mesh data variable 'extra_data' has an " "invalid location='<empty>'."
+            )
+        elif param == "badlocation":
+            match_msg = (
+                "mesh data variable 'extra_data' has an "
+                "invalid location='invalid_location'."
+            )
+        elif param == "baddim":
+            match_msg = (
+                "mesh data variable 'extra_data' does not have the node mesh "
+                "dimension 'node', in its dimensions."
+            )
+        else:
+            raise ValueError(f"unexpected param: {param}")
+
+        with PARSE_UGRID_ON_LOAD.context():
+            with pytest.raises(ValueError, match=match_msg):
+                iris.load(failnc)
diff --git a/lib/iris/tests/unit/fileformats/netcdf/saver/test_Saver__ugrid.py b/lib/iris/tests/unit/fileformats/netcdf/saver/test_Saver__ugrid.py
@@ -508,9 +508,9 @@ def test_multi_cubes_equal_meshes(self):
             self.assertEqual(props["location"], "face")
             self.assertEqual(props["coordinates"], "face_x face_y")
 
-        # the data variables map the appropriate node dimensions
+        # the data variables map the appropriate node dimension
         self.assertEqual(a_props[_VAR_DIMS], ["Mesh2d_faces"])
-        self.assertEqual(b_props[_VAR_DIMS], ["Mesh2d_faces_0"])
+        self.assertEqual(b_props[_VAR_DIMS], ["Mesh2d_faces"])
 
     def test_multi_cubes_different_mesh(self):
         # Check that we can correctly distinguish 2 different meshes.
@@ -1172,7 +1172,50 @@ def test_connectivity_names(self):
             )
             self.assertEqual(expected_names, result_names, fail_msg)
 
-    def _check_two_different_meshes(self, vars):
+    def test_multiple_equal_mesh(self):
+        mesh1 = make_mesh()
+        mesh2 = make_mesh()
+
+        # Save and snapshot the result
+        tempfile_path = self.check_save_mesh([mesh1, mesh2])
+        dims, vars = scan_dataset(tempfile_path)
+
+        # In this case there should be only *one* mesh.
+        mesh_names = vars_meshnames(vars)
+        self.assertEqual(1, len(mesh_names))
+
+        # Check it has the correct number of coords + conns (no duplicates)
+        # Should have 2 each X and Y coords (face+node): _no_ edge coords.
+        coord_vars_x = vars_w_props(vars, standard_name="longitude")
+        coord_vars_y = vars_w_props(vars, standard_name="latitude")
+        self.assertEqual(2, len(coord_vars_x))
+        self.assertEqual(2, len(coord_vars_y))
+
+        # Check the connectivities are all present: _only_ 1 var of each type.
+        for conn in mesh1.all_connectivities:
+            if conn is not None:
+                conn_vars = vars_w_props(vars, cf_role=conn.cf_role)
+                self.assertEqual(1, len(conn_vars))
+
+    def test_multiple_different_meshes(self):
+        # Create 2 meshes with different faces, but same edges.
+        # N.B. they should *not* then share an edge dimension !
+        mesh1 = make_mesh(n_faces=3, n_edges=2)
+        mesh2 = make_mesh(n_faces=4, n_edges=2)
+
+        # Save and snapshot the result
+        tempfile_path = self.check_save_mesh([mesh1, mesh2])
+        dims, vars = scan_dataset(tempfile_path)
+
+        # Check the dims are as expected
+        self.assertEqual(dims["Mesh2d_faces"], 3)
+        self.assertEqual(dims["Mesh2d_faces_0"], 4)
+        # There are no 'second' edge and node dims
+        self.assertEqual(dims["Mesh2d_nodes"], 5)
+        self.assertEqual(dims["Mesh2d_edge"], 2)
+
+        # Check there are two independent meshes in the file...
+
         # there are exactly 2 meshes in the file
         mesh_names = vars_meshnames(vars)
         self.assertEqual(sorted(mesh_names), ["Mesh2d", "Mesh2d_0"])
@@ -1188,15 +1231,15 @@ def _check_two_different_meshes(self, vars):
 
         # mesh2
         self.assertEqual(
-            vars_meshdim(vars, "node", mesh_name="Mesh2d_0"), "Mesh2d_nodes_0"
+            vars_meshdim(vars, "node", mesh_name="Mesh2d_0"), "Mesh2d_nodes"
         )
         self.assertEqual(
             vars_meshdim(vars, "face", mesh_name="Mesh2d_0"), "Mesh2d_faces_0"
         )
         if "edge_coordinates" in vars["Mesh2d_0"]:
             self.assertEqual(
                 vars_meshdim(vars, "edge", mesh_name="Mesh2d_0"),
-                "Mesh2d_edge_0",
+                "Mesh2d_edge",
             )
 
         # the relevant coords + connectivities are also distinct
@@ -1215,63 +1258,19 @@ def _check_two_different_meshes(self, vars):
             )
 
         # mesh2
-        self.assertEqual(vars["node_x_0"][_VAR_DIMS], ["Mesh2d_nodes_0"])
+        self.assertEqual(vars["node_x_0"][_VAR_DIMS], ["Mesh2d_nodes"])
         self.assertEqual(vars["face_x_0"][_VAR_DIMS], ["Mesh2d_faces_0"])
         self.assertEqual(
             vars["mesh2d_faces_0"][_VAR_DIMS],
             ["Mesh2d_faces_0", "Mesh2d_0_face_N_nodes"],
         )
         if "edge_coordinates" in vars["Mesh2d_0"]:
-            self.assertEqual(vars["longitude_0"][_VAR_DIMS], ["Mesh2d_edge_0"])
+            self.assertEqual(vars["longitude_0"][_VAR_DIMS], ["Mesh2d_edge"])
             self.assertEqual(
                 vars["mesh2d_edge_0"][_VAR_DIMS],
-                ["Mesh2d_edge_0", "Mesh2d_0_edge_N_nodes"],
+                ["Mesh2d_edge", "Mesh2d_0_edge_N_nodes"],
             )
 
-    def test_multiple_equal_mesh(self):
-        mesh1 = make_mesh()
-        mesh2 = make_mesh()
-
-        # Save and snapshot the result
-        tempfile_path = self.check_save_mesh([mesh1, mesh2])
-        dims, vars = scan_dataset(tempfile_path)
-
-        # In this case there should be only *one* mesh.
-        mesh_names = vars_meshnames(vars)
-        self.assertEqual(1, len(mesh_names))
-
-        # Check it has the correct number of coords + conns (no duplicates)
-        # Should have 2 each X and Y coords (face+node): _no_ edge coords.
-        coord_vars_x = vars_w_props(vars, standard_name="longitude")
-        coord_vars_y = vars_w_props(vars, standard_name="latitude")
-        self.assertEqual(2, len(coord_vars_x))
-        self.assertEqual(2, len(coord_vars_y))
-
-        # Check the connectivities are all present: _only_ 1 var of each type.
-        for conn in mesh1.all_connectivities:
-            if conn is not None:
-                conn_vars = vars_w_props(vars, cf_role=conn.cf_role)
-                self.assertEqual(1, len(conn_vars))
-
-    def test_multiple_different_meshes(self):
-        # Create 2 meshes with different faces, but same edges.
-        # N.B. they should *not* then share an edge dimension !
-        mesh1 = make_mesh(n_faces=3, n_edges=2)
-        mesh2 = make_mesh(n_faces=4, n_edges=2)
-
-        # Save and snapshot the result
-        tempfile_path = self.check_save_mesh([mesh1, mesh2])
-        dims, vars = scan_dataset(tempfile_path)
-
-        # Check there are two independent meshes
-        self._check_two_different_meshes(vars)
-
-        # Check the dims are as expected
-        self.assertEqual(dims["Mesh2d_faces"], 3)
-        self.assertEqual(dims["Mesh2d_faces_0"], 4)
-        self.assertEqual(dims["Mesh2d_edge"], 2)
-        self.assertEqual(dims["Mesh2d_edge_0"], 2)
-
 
 # WHEN MODIFYING THIS MODULE, CHECK IF ANY CORRESPONDING CHANGES ARE NEEDED IN
 # :mod:`iris.tests.unit.fileformats.netcdf.test_Saver__lazy.`