diff --git a/cf/field.py b/cf/field.py index b571158ea9..0f6b400c68 100644 --- a/cf/field.py +++ b/cf/field.py @@ -6903,7 +6903,8 @@ def collapse( # --------------------------------------------------------- # Update dimension coordinates, auxiliary coordinates, - # cell measures and domain ancillaries + # cell measures, domain ancillaries, domain_topologies, + # and cell connectivities. # --------------------------------------------------------- for axis, domain_axis in collapse_axes.items(): # Ignore axes which are already size 1 @@ -6911,10 +6912,16 @@ def collapse( if size == 1: continue - # REMOVE all cell measures and domain ancillaries - # which span this axis + # REMOVE all cell measures, domain ancillaries, + # domain_topologies, and cell connectivities which + # span this axis c = f.constructs.filter( - filter_by_type=("cell_measure", "domain_ancillary"), + filter_by_type=( + "cell_measure", + "domain_ancillary", + "domain_topology", + "cell_connectivity", + ), filter_by_axis=(axis,), axis_mode="or", todict=True, diff --git a/cf/test/create_test_files.py b/cf/test/create_test_files.py index e2b6cf6f48..604f38c787 100644 --- a/cf/test/create_test_files.py +++ b/cf/test/create_test_files.py @@ -2228,6 +2228,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") @@ -2341,6 +2485,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/cf/test/test_UGRID.py b/cf/test/test_UGRID.py index 7d80ba3166..848da10593 100644 --- a/cf/test/test_UGRID.py +++ b/cf/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 combinations 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 = cf.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 = cf.read(self.filename1) @@ -177,9 +197,103 @@ 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 cf.read and cf.write with UGRID fields.""" + # Face, edge, and point fields that are all part of the same + # UGRID mesh + ugrid = cf.example_fields(8, 9, 10) + face, edge, point = (0, 1, 2) + + tmpfile = "tmpfileu.nc" + # Test for equality with the fields defined in memory. Only + # works for face and edge fields. + for cell in (face, edge): + f = ugrid[cell] + cf.write(f, tmpfile) + g = cf.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]) + + cf.write(f, tmpfile) + + # Check that there's only one mesh variable in the file + self.assertEqual(n_mesh_variables(tmpfile), 1) + + g = cf.read(tmpfile) + self.assertEqual(len(g), len(f)) + + cf.write(g, tmpfile1) + + # Check that there's only one mesh variable in the file + self.assertEqual(n_mesh_variables(tmpfile1), 1) + + h = cf.read(tmpfile1) + self.assertEqual(len(h), len(g)) + self.assertTrue(h[0].equals(g[0])) + + def test_read_write_UGRID_domain(self): + """Test the cf.read and cf.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 cf.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] + cf.write(d, tmpfile) + e = cf.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]) + + cf.write(d, tmpfile) + + # Check that there's only one mesh variable in the file + self.assertEqual(n_mesh_variables(tmpfile), 1) + + e = cf.read(tmpfile, domain=True) + + self.assertEqual(len(e), len(d)) + + cf.write(e, tmpfile1) + + # Check that there's only one mesh variable in the file + self.assertEqual(n_mesh_variables(tmpfile1), 1) + + f = cf.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 = cf.read(self.filename3, domain=True) + face, edge, point = (2, 1, 0) if __name__ == "__main__": diff --git a/cf/test/test_collapse.py b/cf/test/test_collapse.py index c7c866c8b6..0f987fb969 100644 --- a/cf/test/test_collapse.py +++ b/cf/test/test_collapse.py @@ -1,6 +1,8 @@ +import atexit import datetime import faulthandler import os +import tempfile import unittest import numpy @@ -9,6 +11,25 @@ import cf +n_tmpfiles = 1 +tmpfiles = [ + tempfile.mkstemp("_test_collapse.nc", dir=os.getcwd())[1] + for i in range(n_tmpfiles) +] +[tmpfile] = tmpfiles + + +def _remove_tmpfiles(): + """Try to remove defined temporary files by deleting their paths.""" + for f in tmpfiles: + try: + os.remove(f) + except OSError: + pass + + +atexit.register(_remove_tmpfiles) + class Field_collapseTest(unittest.TestCase): def setUp(self): @@ -791,6 +812,19 @@ def test_Field_collapse_non_positive_weights(self): # compute time g.array + def test_Field_collapse_ugrid(self): + """Check that UGRID constructs are removed after collapsing.""" + f = cf.example_field(8) + self.assertTrue(f.domain_topologies()) + self.assertTrue(f.cell_connectivities()) + + f = f.collapse("area: mean") + self.assertFalse(f.domain_topologies()) + self.assertFalse(f.cell_connectivities()) + + # Check the collpsed fields writes + cf.write(f, tmpfile) + if __name__ == "__main__": print("Run date:", datetime.datetime.now()) diff --git a/cf/test/test_read_write.py b/cf/test/test_read_write.py index 328be3e3c7..56a554d771 100644 --- a/cf/test/test_read_write.py +++ b/cf/test/test_read_write.py @@ -346,9 +346,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 @@ -426,9 +425,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. cf.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 = cf.example_fields(0, 1, 2, 3, 4, 5, 6, 7) del append_ex_fields[1] # note: can remove after Issue #141 closed if fmt in "NETCDF4_CLASSIC":