Mesh full comparison (#4439)

* Fixes to transposed Connectivity equality.

* Fixes to transposed Connectivity equality.

* Make _DimensionalMetadata and Connectivity equality preserve laziness.

* Full Mesh equality check.

* Remove outdated comment about Mesh equality.

* What's new entry.

* Extra What's New.

* iris.util.array_equal fixes

* Dask optimise array_equal.

* Add a metadata difference to test_Mesh test_different.

* Adapt tests to improved Mesh equality.

* What's New correction.

* Give Stephen UGRID credit too.

Co-authored-by: Bill Little <bill.james.little@gmail.com>

docs/src/whatsnew/latest.rst

@@ -34,11 +34,11 @@ This document explains the changes made to Iris for this release
 ✨ Features
 ===========
 
-#. `@bjlittle`_, `@pp-mo`_ and `@trexfeathers`_ added support for unstructured
-   meshes, as described by `UGRID`_. This involved adding a data model (:pull:`3968`,
-   :pull:`4014`, :pull:`4027`, :pull:`4036`, :pull:`4053`) and API (:pull:`4063`,
-   :pull:`4064`), and supporting representation (:pull:`4033`, :pull:`4054`) of
-   data on meshes.
+#. `@bjlittle`_, `@pp-mo`_, `@trexfeathers`_ and `@stephenworsley`_ added
+   support for unstructured meshes, as described by `UGRID`_. This involved
+   adding a data model (:pull:`3968`, :pull:`4014`, :pull:`4027`, :pull:`4036`,
+   :pull:`4053`, :pull:`4439`) and API (:pull:`4063`, :pull:`4064`), and
+   supporting representation (:pull:`4033`, :pull:`4054`) of data on meshes.
    Most of this new API can be found in :mod:`iris.experimental.ugrid`. The key
    objects introduced are :class:`iris.experimental.ugrid.mesh.Mesh`,
    :class:`iris.experimental.ugrid.mesh.MeshCoord` and
@@ -130,6 +130,11 @@ This document explains the changes made to Iris for this release
 #. `@wjbenfold`_ changed :meth:`iris.util.points_step` to stop it from warning
    when applied to a single point (:issue:`4250`, :pull:`4367`)
 
+#. `@trexfeathers`_ changed :class:`~iris.coords._DimensionalMetadata` and
+   :class:`~iris.experimental.ugrid.Connectivity` equality methods to preserve
+   array laziness, allowing efficient comparisons even with larger-than-memory
+   objects. (:pull:`4439`)
+
 
 💣 Incompatible Changes
 =======================

lib/iris/coords.py

@@ -343,27 +343,27 @@ def __repr__(self):
 
     def __eq__(self, other):
         # Note: this method includes bounds handling code, but it only runs
-        # within Coord type instances, as only these allow bounds to be set.
+        #  within Coord type instances, as only these allow bounds to be set.
 
         eq = NotImplemented
         # If the other object has a means of getting its definition, then do
-        # the  comparison, otherwise return a NotImplemented to let Python try
-        # to resolve the operator elsewhere.
+        #  the comparison, otherwise return a NotImplemented to let Python try
+        #  to resolve the operator elsewhere.
         if hasattr(other, "metadata"):
             # metadata comparison
             eq = self.metadata == other.metadata
             # data values comparison
             if eq and eq is not NotImplemented:
                 eq = iris.util.array_equal(
-                    self._values, other._values, withnans=True
+                    self._core_values(), other._core_values(), withnans=True
                 )
 
             # Also consider bounds, if we have them.
             # (N.B. though only Coords can ever actually *have* bounds).
             if eq and eq is not NotImplemented:
                 if self.has_bounds() and other.has_bounds():
                     eq = iris.util.array_equal(
-                        self.bounds, other.bounds, withnans=True
+                        self.core_bounds(), other.core_bounds(), withnans=True
                     )
                 else:
                     eq = not self.has_bounds() and not other.has_bounds()

lib/iris/cube.py

@@ -3507,6 +3507,8 @@ def __eq__(self, other):
 
             # Having checked everything else, check approximate data equality.
             if result:
+                # TODO: why do we use allclose() here, but strict equality in
+                #  _DimensionalMetadata (via util.array_equal())?
                 result = da.allclose(
                     self.core_data(), other.core_data()
                 ).compute()

lib/iris/experimental/ugrid/mesh.py

@@ -27,7 +27,7 @@
 from ...config import get_logger
 from ...coords import AuxCoord, _DimensionalMetadata
 from ...exceptions import ConnectivityNotFoundError, CoordinateNotFoundError
-from ...util import guess_coord_axis
+from ...util import array_equal, guess_coord_axis
 from .metadata import ConnectivityMetadata, MeshCoordMetadata, MeshMetadata
 
 # Configure the logger.
@@ -483,12 +483,19 @@ 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()
-                    ).all()
+                        self.shape == other.shape
+                        and self.src_dim == other.src_dim
+                    ) or (
+                        self.shape == other.shape[::-1]
+                        and self.src_dim == other.tgt_dim
+                    )
+                if eq:
+                    eq = array_equal(
+                        self.indices_by_src(self.core_indices()),
+                        other.indices_by_src(other.core_indices()),
+                    )
         return eq
 
     def transpose(self):
@@ -939,8 +946,16 @@ def axes_assign(coord_list):
         return cls(**mesh_kwargs)
 
     def __eq__(self, other):
-        # TBD: this is a minimalist implementation and requires to be revisited
-        return id(self) == id(other)
+        result = NotImplemented
+
+        if isinstance(other, Mesh):
+            result = self.metadata == other.metadata
+            if result:
+                result = self.all_coords == other.all_coords
+            if result:
+                result = self.all_connectivities == other.all_connectivities
+
+        return result
 
     def __hash__(self):
         # Allow use in sets and as dictionary keys, as is done for :class:`iris.cube.Cube`.
@@ -2883,9 +2898,7 @@ def copy(self, points=None, bounds=None):
         """
         # Override Coord.copy, so that we can ensure it does not duplicate the
         # Mesh object (via deepcopy).
-        # This avoids copying Meshes.  It is also required to allow a copied
-        # MeshCoord to be == the original, since for now Mesh == is only true
-        # for the same identical object.
+        # This avoids copying Meshes.
 
         # FOR NOW: also disallow changing points/bounds at all.
         if points is not None or bounds is not None:

lib/iris/tests/unit/cube/test_Cube.py

@@ -2351,13 +2351,19 @@ def test_fail_meshcoords_different_locations(self):
                 aux_coords_and_dims=[(meshco_1, 0), (meshco_2, 0)],
             )
 
+    def test_meshcoords_equal_meshes(self):
+        meshco_x = sample_meshcoord(axis="x")
+        meshco_y = sample_meshcoord(axis="y")
+        n_faces = meshco_x.shape[0]
+        Cube(
+            np.zeros(n_faces),
+            aux_coords_and_dims=[(meshco_x, 0), (meshco_y, 0)],
+        )
+
     def test_fail_meshcoords_different_meshes(self):
-        # Same as successful 'multi_mesh', but not sharing the same mesh.
-        # This one *is* an error.
-        # But that could relax in future, if we allow mesh equality testing
-        # (i.e. "mesh_a == mesh_b" when not "mesh_a is mesh_b")
         meshco_x = sample_meshcoord(axis="x")
         meshco_y = sample_meshcoord(axis="y")  # Own (different) mesh
+        meshco_y.mesh.long_name = "new_name"
         n_faces = meshco_x.shape[0]
         with self.assertRaisesRegex(ValueError, "Mesh.* does not match"):
             Cube(
@@ -2412,11 +2418,20 @@ def test_add_multiple(self):
         cube.add_aux_coord(new_meshco_y, 1)
         self.assertEqual(len(cube.coords(mesh_coords=True)), 3)
 
+    def test_add_equal_mesh(self):
+        # Make a duplicate y-meshco, and rename so it can add into the cube.
+        cube = self.cube
+        # Create 'meshco_y' duplicate, but a new mesh
+        meshco_y = sample_meshcoord(axis="y")
+        cube.add_aux_coord(meshco_y, 1)
+        self.assertIn(meshco_y, cube.coords(mesh_coords=True))
+
     def test_fail_different_mesh(self):
         # Make a duplicate y-meshco, and rename so it can add into the cube.
         cube = self.cube
         # Create 'meshco_y' duplicate, but a new mesh
         meshco_y = sample_meshcoord(axis="y")
+        meshco_y.mesh.long_name = "new_name"
         msg = "does not match existing cube mesh"
         with self.assertRaisesRegex(ValueError, msg):
             cube.add_aux_coord(meshco_y, 1)
@@ -2481,17 +2496,18 @@ def test_copied_cube_match(self):
         cube2 = cube.copy()
         self.assertEqual(cube, cube2)
 
-    def test_same_mesh_match(self):
+    def test_equal_mesh_match(self):
         cube1 = self.cube
         # re-create an identical cube, using the same mesh.
-        _add_test_meshcube(self, mesh=self.mesh)
+        _add_test_meshcube(self)
         cube2 = self.cube
         self.assertEqual(cube1, cube2)
 
     def test_new_mesh_different(self):
         cube1 = self.cube
-        # re-create an identical cube, using the same mesh.
+        # re-create an identical cube, using a different mesh.
         _add_test_meshcube(self)
+        self.cube.mesh.long_name = "new_name"
         cube2 = self.cube
         self.assertNotEqual(cube1, cube2)
 

lib/iris/tests/unit/experimental/ugrid/mesh/test_Connectivity.py

@@ -23,7 +23,7 @@ def setUp(self):
         # Crete an instance, with non-default arguments to allow testing of
         # correct property setting.
         self.kwargs = {
-            "indices": np.linspace(1, 9, 9, dtype=int).reshape((3, -1)),
+            "indices": np.linspace(1, 12, 12, dtype=int).reshape((4, -1)),
             "cf_role": "face_node_connectivity",
             "long_name": "my_face_nodes",
             "var_name": "face_nodes",
@@ -91,7 +91,7 @@ def test_lazy_src_lengths(self):
         self.assertTrue(is_lazy_data(self.connectivity.lazy_src_lengths()))
 
     def test_src_lengths(self):
-        expected = [3, 3, 3]
+        expected = [4, 4, 4]
         self.assertArrayEqual(expected, self.connectivity.src_lengths())
 
     def test___str__(self):
@@ -102,7 +102,7 @@ def test___str__(self):
 
     def test___repr__(self):
         expected = (
-            "Connectivity(array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]), "
+            "Connectivity(array([[ 1,  2,  3], [ 4,  5,  6], [ 7,  8,  9], [10, 11, 12]]), "
             "cf_role='face_node_connectivity', long_name='my_face_nodes', "
             "var_name='face_nodes', attributes={'notes': 'this is a test'}, "
             "start_index=1, src_dim=1)"
@@ -122,7 +122,7 @@ def test___eq__(self):
         equivalent_kwargs["src_dim"] = 1 - self.kwargs["src_dim"]
         equivalent = Connectivity(**equivalent_kwargs)
         self.assertFalse(
-            (equivalent.indices == self.connectivity.indices).all()
+            np.array_equal(equivalent.indices, self.connectivity.indices)
         )
         self.assertEqual(equivalent, self.connectivity)
 

lib/iris/tests/unit/experimental/ugrid/mesh/test_Mesh.py

@@ -81,7 +81,7 @@ def setUpClass(cls):
         cls.kwargs = {
             "topology_dimension": 1,
             "node_coords_and_axes": ((cls.NODE_LON, "x"), (cls.NODE_LAT, "y")),
-            "connectivities": cls.EDGE_NODE,
+            "connectivities": [cls.EDGE_NODE],
             "long_name": "my_topology_mesh",
             "var_name": "mesh",
             "attributes": {"notes": "this is a test"},
@@ -124,6 +124,34 @@ def test___repr__(self):
         )
         self.assertEqual(expected, self.mesh.__repr__())
 
+    def test___eq__(self):
+        # The dimension names do not participate in equality.
+        equivalent_kwargs = self.kwargs.copy()
+        equivalent_kwargs["node_dimension"] = "something_else"
+        equivalent = mesh.Mesh(**equivalent_kwargs)
+        self.assertEqual(equivalent, self.mesh)
+
+    def test_different(self):
+        different_kwargs = self.kwargs.copy()
+        different_kwargs["long_name"] = "new_name"
+        different = mesh.Mesh(**different_kwargs)
+        self.assertNotEqual(different, self.mesh)
+
+        different_kwargs = self.kwargs.copy()
+        ncaa = self.kwargs["node_coords_and_axes"]
+        new_lat = ncaa[1][0].copy(points=ncaa[1][0].points + 1)
+        new_ncaa = (ncaa[0], (new_lat, "y"))
+        different_kwargs["node_coords_and_axes"] = new_ncaa
+        different = mesh.Mesh(**different_kwargs)
+        self.assertNotEqual(different, self.mesh)
+
+        different_kwargs = self.kwargs.copy()
+        conns = self.kwargs["connectivities"]
+        new_conn = conns[0].copy(conns[0].indices + 1)
+        different_kwargs["connectivities"] = new_conn
+        different = mesh.Mesh(**different_kwargs)
+        self.assertNotEqual(different, self.mesh)
+
     def test_all_connectivities(self):
         expected = mesh.Mesh1DConnectivities(self.EDGE_NODE)
         self.assertEqual(expected, self.mesh.all_connectivities)

lib/iris/tests/unit/experimental/ugrid/mesh/test_MeshCoord.py

@@ -183,19 +183,24 @@ def setUp(self):
     def _create_common_mesh(self, **kwargs):
         return sample_meshcoord(mesh=self.mesh, **kwargs)
 
-    def test_same_mesh(self):
+    def test_identical_mesh(self):
         meshcoord1 = self._create_common_mesh()
         meshcoord2 = self._create_common_mesh()
         self.assertEqual(meshcoord2, meshcoord1)
 
-    def test_different_identical_mesh(self):
-        # For equality, must have the SAME mesh (at present).
+    def test_equal_mesh(self):
         mesh1 = sample_mesh()
-        mesh2 = sample_mesh()  # Presumably identical, but not the same
+        mesh2 = sample_mesh()
+        meshcoord1 = sample_meshcoord(mesh=mesh1)
+        meshcoord2 = sample_meshcoord(mesh=mesh2)
+        self.assertEqual(meshcoord2, meshcoord1)
+
+    def test_different_mesh(self):
+        mesh1 = sample_mesh()
+        mesh2 = sample_mesh()
+        mesh2.long_name = "new_name"
         meshcoord1 = sample_meshcoord(mesh=mesh1)
         meshcoord2 = sample_meshcoord(mesh=mesh2)
-        # These should NOT compare, because the Meshes are not identical : at
-        # present, Mesh equality is not implemented (i.e. limited to identity)
         self.assertNotEqual(meshcoord2, meshcoord1)
 
     def test_different_location(self):

lib/iris/tests/unit/fileformats/netcdf/test_Saver__ugrid.py

@@ -487,7 +487,7 @@ def test_multi_cubes_different_locations(self):
         self.assertEqual(v_a[_VAR_DIMS], [face_dim])
         self.assertEqual(v_b[_VAR_DIMS], [node_dim])
 
-    def test_multi_cubes_identical_meshes(self):
+    def test_multi_cubes_equal_meshes(self):
         # Make 2 identical meshes
         # NOTE: *can't* name these explicitly, as it stops them being identical.
         mesh1 = make_mesh()
@@ -499,34 +499,27 @@ def test_multi_cubes_identical_meshes(self):
         tempfile_path = self.check_save_cubes([cube1, cube2])
         dims, vars = scan_dataset(tempfile_path)
 
-        # there are exactly 2 meshes in the file
+        # there is exactly 1 mesh in the file
         mesh_names = vars_meshnames(vars)
-        self.assertEqual(sorted(mesh_names), ["Mesh2d", "Mesh2d_0"])
+        self.assertEqual(sorted(mesh_names), ["Mesh2d"])
 
-        # they use different dimensions
+        # same dimensions
         self.assertEqual(
             vars_meshdim(vars, "node", mesh_name="Mesh2d"), "Mesh2d_nodes"
         )
         self.assertEqual(
             vars_meshdim(vars, "face", mesh_name="Mesh2d"), "Mesh2d_faces"
         )
-        self.assertEqual(
-            vars_meshdim(vars, "node", mesh_name="Mesh2d_0"), "Mesh2d_nodes_0"
-        )
-        self.assertEqual(
-            vars_meshdim(vars, "face", mesh_name="Mesh2d_0"), "Mesh2d_faces_0"
-        )
 
         # there are exactly two data-variables with a 'mesh' property
         mesh_datavars = vars_w_props(vars, mesh="*")
         self.assertEqual(["a", "b"], list(mesh_datavars))
 
-        # the data variables reference the two separate meshes
+        # the data variables reference the same mesh
         a_props, b_props = vars["a"], vars["b"]
-        self.assertEqual(a_props["mesh"], "Mesh2d")
-        self.assertEqual(a_props["location"], "face")
-        self.assertEqual(b_props["mesh"], "Mesh2d_0")
-        self.assertEqual(b_props["location"], "face")
+        for props in a_props, b_props:
+            self.assertEqual(props["mesh"], "Mesh2d")
+            self.assertEqual(props["location"], "face")
 
         # the data variables map the appropriate node dimensions
         self.assertEqual(a_props[_VAR_DIMS], ["Mesh2d_faces"])
@@ -1234,24 +1227,14 @@ def _check_two_different_meshes(self, vars):
                 ["Mesh2d_edge_0", "Mesh2d_0_edge_N_nodes"],
             )
 
-    def test_multiple_identical_meshes(self):
+    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)
 
-        # Check there are two independent meshes
-        self._check_two_different_meshes(vars)
-
-    def test_multiple_same_mesh(self):
-        mesh = make_mesh()
-
-        # Save and snapshot the result
-        tempfile_path = self.check_save_mesh([mesh, mesh])
-        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))
@@ -1264,7 +1247,7 @@ def test_multiple_same_mesh(self):
         self.assertEqual(2, len(coord_vars_y))
 
         # Check the connectivities are all present: _only_ 1 var of each type.
-        for conn in mesh.all_connectivities:
+        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))