Mesh saveload fix

docs/src/whatsnew/latest.rst

@@ -43,6 +43,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
 =======================
 

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

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":

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):

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)

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.`