From a6d80b817dc37ebacfc8060a4d962c14ba8f60ff Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 20 Jun 2024 17:19:32 +0100 Subject: [PATCH 1/5] MeshCoord coord-system comes from location coords. Fix MeshCoord.copy() to respect that. --- lib/iris/experimental/ugrid/mesh.py | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/lib/iris/experimental/ugrid/mesh.py b/lib/iris/experimental/ugrid/mesh.py index a798f7af77..d8a62e047c 100644 --- a/lib/iris/experimental/ugrid/mesh.py +++ b/lib/iris/experimental/ugrid/mesh.py @@ -2781,6 +2781,11 @@ def fix_repr(val): ) raise ValueError(msg) + # Don't use 'coord_system' as a constructor arg, since for + # MeshCoords it is deduced from the mesh. + # (Otherwise a non-None coord_system breaks the 'copy' operation) + use_metadict.pop("coord_system") + # Call parent constructor to handle the common constructor args. super().__init__(points, bounds=bounds, **use_metadict) @@ -2806,8 +2811,19 @@ def axis(self): @property def coord_system(self): - """The coordinate-system of a MeshCoord is always 'None'.""" - return None + """The coordinate-system of a MeshCoord comes from its origin coords.""" + select_kwargs = { + f"include_{self.location}s": True, + "axis": self.axis, + } + try: + location_coord = self.mesh.coord(**select_kwargs) + coord_system = location_coord.coord_system + except CoordinateNotFoundError: + # No location coord : possible in UGRID but probably not Iris (at present). + coord_system = None + + return coord_system @coord_system.setter def coord_system(self, value): From 0b2c4abe05cfedb55dd5ac0123a5affc217915d7 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Wed, 26 Jun 2024 13:44:00 +0100 Subject: [PATCH 2/5] Add extra context in comments. --- lib/iris/experimental/ugrid/mesh.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/lib/iris/experimental/ugrid/mesh.py b/lib/iris/experimental/ugrid/mesh.py index d8a62e047c..70eb20aa15 100644 --- a/lib/iris/experimental/ugrid/mesh.py +++ b/lib/iris/experimental/ugrid/mesh.py @@ -2811,16 +2811,22 @@ def axis(self): @property def coord_system(self): - """The coordinate-system of a MeshCoord comes from its origin coords.""" + """The coordinate-system of a MeshCoord comes from the related mesh coord.""" + # This matches where the coord metadata is drawn from. + # See : https://github.com/SciTools/iris/issues/4860 select_kwargs = { f"include_{self.location}s": True, "axis": self.axis, } try: + # NOTE: at present, a MeshCoord *always* references a related mesh coord of + # its location, from which it's points are taken. + # However this might change in future .. + # see : https://github.com/SciTools/iris/discussions/4438#bounds-no-points location_coord = self.mesh.coord(**select_kwargs) coord_system = location_coord.coord_system except CoordinateNotFoundError: - # No location coord : possible in UGRID but probably not Iris (at present). + # No such coord : possible in UGRID, but probably not Iris (at present). coord_system = None return coord_system From e418058464d6ef8e9cafc6393d80561d0c82ed0e Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Mon, 15 Jul 2024 18:42:52 +0100 Subject: [PATCH 3/5] Test meshcoord coord-system derivation. --- .../experimental/test_meshcoord_coordsys.py | 112 ++++++++++++++++++ 1 file changed, 112 insertions(+) create mode 100644 lib/iris/tests/integration/experimental/test_meshcoord_coordsys.py diff --git a/lib/iris/tests/integration/experimental/test_meshcoord_coordsys.py b/lib/iris/tests/integration/experimental/test_meshcoord_coordsys.py new file mode 100644 index 0000000000..30500e8b02 --- /dev/null +++ b/lib/iris/tests/integration/experimental/test_meshcoord_coordsys.py @@ -0,0 +1,112 @@ +# 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. +"""Integration tests for MeshCoord.coord_system() behaviour.""" + +import pytest + +import iris +from iris.coord_systems import GeogCS +from iris.experimental.ugrid.load import PARSE_UGRID_ON_LOAD +from iris.tests.stock.netcdf import ncgen_from_cdl + +TEST_CDL = """ +netcdf mesh_test { + dimensions: + node = 3 ; + face = 1 ; + vertex = 3 ; + 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 ; + float node_data(node) ; + node_data:coordinates = "node_x node_y" ; + node_data:location = "node" ; + node_data:mesh = "mesh" ; +} +""" + +BASIC_CRS_STR = """ + int crs ; + crs:grid_mapping_name = "latitude_longitude" ; + crs:semi_major_axis = 6371000.0 ; +""" + + +def find_i_line(lines, match): + (i_line,) = [i for i, line in enumerate(lines) if match in line] + return i_line + + +def make_file(nc_path, node_x_crs=False, node_y_crs=False): + lines = TEST_CDL.split("\n") + + includes = ["node_x"] if node_x_crs else [] + includes += ["node_y"] if node_y_crs else [] + includes = " ".join(includes) + if includes: + # insert a grid-mapping + i_line = find_i_line(lines, "variables:") + lines[i_line + 1 : i_line + 1] = BASIC_CRS_STR.split("\n") + + i_line = find_i_line(lines, "float node_data(") + ref_lines = [ + # NOTE: space to match the surrounding indent + f' node_data:grid_mapping = "crs: {includes}" ;', + # NOTE: this is already present + # f'node_data:coordinates = "{includes}" ;' + ] + lines[i_line + 1 : i_line + 1] = ref_lines + + cdl_str = "\n".join(lines) + ncgen_from_cdl(cdl_str=cdl_str, cdl_path=None, nc_path=nc_path) + + +@pytest.mark.parametrize("cs_axes", ["--", "x-", "-y", "xy"]) +def test_default_mesh_cs(tmp_path, cs_axes): + """Test coord-systems of mesh cube and coords, if location coords have a crs.""" + nc_path = tmp_path / "test_temp.nc" + do_x = "x" in cs_axes + do_y = "y" in cs_axes + make_file(nc_path, node_x_crs=do_x, node_y_crs=do_y) + with PARSE_UGRID_ON_LOAD.context(): + cube = iris.load_cube(nc_path, "node_data") + meshco_x, meshco_y = [cube.coord(mesh_coords=True, axis=ax) for ax in ("x", "y")] + # NOTE: at present, none of these load with a coordinate system, + # because we don't support the extended grid-mapping syntax. + # see: https://github.com/SciTools/iris/issues/3388 + assert meshco_x.coord_system is None + assert meshco_y.coord_system is None + + +def test_assigned_mesh_cs(tmp_path): + # Check that when a coord system is manually assigned to a location coord, + # the corresponding meshcoord reports the same cs. + nc_path = tmp_path / "test_temp.nc" + make_file(nc_path) + with PARSE_UGRID_ON_LOAD.context(): + cube = iris.load_cube(nc_path, "node_data") + nodeco_x = cube.mesh.coord(include_nodes=True, axis="x") + meshco_x, meshco_y = [cube.coord(axis=ax) for ax in ("x", "y")] + assert nodeco_x.coord_system is None + assert meshco_x.coord_system is None + assert meshco_y.coord_system is None + assigned_cs = GeogCS(1.0) + nodeco_x.coord_system = assigned_cs + assert meshco_x.coord_system is assigned_cs + assert meshco_y.coord_system is None + # This also affects cube.coord_system(), even though it is an auxcoord, + # since there are no dim-coords, or any other coord with a c-s. + # TODO: this may be a mistake -- see https://github.com/SciTools/iris/issues/6051 + assert cube.coord_system() is assigned_cs From 5cc17b69d60c982bcf44cabfd21af3413fb7e194 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Mon, 15 Jul 2024 18:53:26 +0100 Subject: [PATCH 4/5] Added whatsnew. --- docs/src/whatsnew/latest.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/src/whatsnew/latest.rst b/docs/src/whatsnew/latest.rst index a6915c4170..76920c5ad5 100644 --- a/docs/src/whatsnew/latest.rst +++ b/docs/src/whatsnew/latest.rst @@ -39,6 +39,12 @@ This document explains the changes made to Iris for this release the :class:`~iris.cube.Cube` :attr:`~iris.cube.Cube.mesh_dim` (see :ref:`cube-statistics-collapsing`). (:issue:`5377`, :pull:`6003`) +#. `@pp-mo`_ made a MeshCoord inherit a coordinate system from its location coord, + as it does its metadata. N.B. mesh location coords can not however load a + coordinate system from netcdf at present, as this needs the 'extended' + grid-mappping syntax -- see : :issue:`3388`. + (:issue:`5562`, :pull:`6016`) + 🐛 Bugs Fixed ============= From f8ef8e0bcd4dd8139d4adfc18f00a43653f9a6bf Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Wed, 17 Jul 2024 13:43:03 +0100 Subject: [PATCH 5/5] Review changes. --- lib/iris/experimental/ugrid/mesh.py | 9 ++++++--- .../experimental/test_meshcoord_coordsys.py | 17 +++++++++++++++++ 2 files changed, 23 insertions(+), 3 deletions(-) diff --git a/lib/iris/experimental/ugrid/mesh.py b/lib/iris/experimental/ugrid/mesh.py index c1d77c34ba..651906b748 100644 --- a/lib/iris/experimental/ugrid/mesh.py +++ b/lib/iris/experimental/ugrid/mesh.py @@ -2814,7 +2814,10 @@ def axis(self): @property def coord_system(self): - """The coordinate-system of a MeshCoord comes from the related mesh coord.""" + """The coordinate-system of a MeshCoord. + + It comes from the `related` location coordinate in the mesh. + """ # This matches where the coord metadata is drawn from. # See : https://github.com/SciTools/iris/issues/4860 select_kwargs = { @@ -2822,8 +2825,8 @@ def coord_system(self): "axis": self.axis, } try: - # NOTE: at present, a MeshCoord *always* references a related mesh coord of - # its location, from which it's points are taken. + # NOTE: at present, a MeshCoord *always* references the relevant location + # coordinate in the mesh, from which its points are taken. # However this might change in future .. # see : https://github.com/SciTools/iris/discussions/4438#bounds-no-points location_coord = self.mesh.coord(**select_kwargs) diff --git a/lib/iris/tests/integration/experimental/test_meshcoord_coordsys.py b/lib/iris/tests/integration/experimental/test_meshcoord_coordsys.py index 30500e8b02..353d56004f 100644 --- a/lib/iris/tests/integration/experimental/test_meshcoord_coordsys.py +++ b/lib/iris/tests/integration/experimental/test_meshcoord_coordsys.py @@ -110,3 +110,20 @@ def test_assigned_mesh_cs(tmp_path): # since there are no dim-coords, or any other coord with a c-s. # TODO: this may be a mistake -- see https://github.com/SciTools/iris/issues/6051 assert cube.coord_system() is assigned_cs + + +def test_meshcoord_coordsys_copy(tmp_path): + # Check that copying a meshcoord with a coord system works properly. + nc_path = tmp_path / "test_temp.nc" + make_file(nc_path) + with PARSE_UGRID_ON_LOAD.context(): + cube = iris.load_cube(nc_path, "node_data") + node_coord = cube.mesh.coord(include_nodes=True, axis="x") + assigned_cs = GeogCS(1.0) + node_coord.coord_system = assigned_cs + mesh_coord = cube.coord(axis="x") + assert mesh_coord.coord_system is assigned_cs + meshco_copy = mesh_coord.copy() + assert meshco_copy == mesh_coord + # Note: still the same object, because it is derived from the same node_coord + assert meshco_copy.coord_system is assigned_cs