Skip to content

Commit 17dd34a

Browse files
Allow MeshCoord to have a coord-system (#6016)
* MeshCoord coord-system comes from location coords. Fix MeshCoord.copy() to respect that. * Add extra context in comments. * Test meshcoord coord-system derivation. * Added whatsnew. * Review changes. --------- Co-authored-by: stephenworsley <[email protected]>
1 parent 8604bb8 commit 17dd34a

File tree

3 files changed

+162
-2
lines changed

3 files changed

+162
-2
lines changed

docs/src/whatsnew/latest.rst

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,12 @@ This document explains the changes made to Iris for this release
3939
the :class:`~iris.cube.Cube` :attr:`~iris.cube.Cube.mesh_dim` (see
4040
:ref:`cube-statistics-collapsing`). (:issue:`5377`, :pull:`6003`)
4141

42+
#. `@pp-mo`_ made a MeshCoord inherit a coordinate system from its location coord,
43+
as it does its metadata. N.B. mesh location coords can not however load a
44+
coordinate system from netcdf at present, as this needs the 'extended'
45+
grid-mappping syntax -- see : :issue:`3388`.
46+
(:issue:`5562`, :pull:`6016`)
47+
4248

4349
🐛 Bugs Fixed
4450
=============

lib/iris/experimental/ugrid/mesh.py

Lines changed: 27 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2784,6 +2784,11 @@ def fix_repr(val):
27842784
)
27852785
raise ValueError(msg)
27862786

2787+
# Don't use 'coord_system' as a constructor arg, since for
2788+
# MeshCoords it is deduced from the mesh.
2789+
# (Otherwise a non-None coord_system breaks the 'copy' operation)
2790+
use_metadict.pop("coord_system")
2791+
27872792
# Call parent constructor to handle the common constructor args.
27882793
super().__init__(points, bounds=bounds, **use_metadict)
27892794

@@ -2809,8 +2814,28 @@ def axis(self):
28092814

28102815
@property
28112816
def coord_system(self):
2812-
"""The coordinate-system of a MeshCoord is always 'None'."""
2813-
return None
2817+
"""The coordinate-system of a MeshCoord.
2818+
2819+
It comes from the `related` location coordinate in the mesh.
2820+
"""
2821+
# This matches where the coord metadata is drawn from.
2822+
# See : https://github.com/SciTools/iris/issues/4860
2823+
select_kwargs = {
2824+
f"include_{self.location}s": True,
2825+
"axis": self.axis,
2826+
}
2827+
try:
2828+
# NOTE: at present, a MeshCoord *always* references the relevant location
2829+
# coordinate in the mesh, from which its points are taken.
2830+
# However this might change in future ..
2831+
# see : https://github.com/SciTools/iris/discussions/4438#bounds-no-points
2832+
location_coord = self.mesh.coord(**select_kwargs)
2833+
coord_system = location_coord.coord_system
2834+
except CoordinateNotFoundError:
2835+
# No such coord : possible in UGRID, but probably not Iris (at present).
2836+
coord_system = None
2837+
2838+
return coord_system
28142839

28152840
@coord_system.setter
28162841
def coord_system(self, value):
Lines changed: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
1+
# Copyright Iris contributors
2+
#
3+
# This file is part of Iris and is released under the BSD license.
4+
# See LICENSE in the root of the repository for full licensing details.
5+
"""Integration tests for MeshCoord.coord_system() behaviour."""
6+
7+
import pytest
8+
9+
import iris
10+
from iris.coord_systems import GeogCS
11+
from iris.experimental.ugrid.load import PARSE_UGRID_ON_LOAD
12+
from iris.tests.stock.netcdf import ncgen_from_cdl
13+
14+
TEST_CDL = """
15+
netcdf mesh_test {
16+
dimensions:
17+
node = 3 ;
18+
face = 1 ;
19+
vertex = 3 ;
20+
variables:
21+
int mesh ;
22+
mesh:cf_role = "mesh_topology" ;
23+
mesh:topology_dimension = 2 ;
24+
mesh:node_coordinates = "node_x node_y" ;
25+
mesh:face_node_connectivity = "face_nodes" ;
26+
float node_x(node) ;
27+
node_x:standard_name = "longitude" ;
28+
float node_y(node) ;
29+
node_y:standard_name = "latitude" ;
30+
int face_nodes(face, vertex) ;
31+
face_nodes:cf_role = "face_node_connectivity" ;
32+
face_nodes:start_index = 0 ;
33+
float node_data(node) ;
34+
node_data:coordinates = "node_x node_y" ;
35+
node_data:location = "node" ;
36+
node_data:mesh = "mesh" ;
37+
}
38+
"""
39+
40+
BASIC_CRS_STR = """
41+
int crs ;
42+
crs:grid_mapping_name = "latitude_longitude" ;
43+
crs:semi_major_axis = 6371000.0 ;
44+
"""
45+
46+
47+
def find_i_line(lines, match):
48+
(i_line,) = [i for i, line in enumerate(lines) if match in line]
49+
return i_line
50+
51+
52+
def make_file(nc_path, node_x_crs=False, node_y_crs=False):
53+
lines = TEST_CDL.split("\n")
54+
55+
includes = ["node_x"] if node_x_crs else []
56+
includes += ["node_y"] if node_y_crs else []
57+
includes = " ".join(includes)
58+
if includes:
59+
# insert a grid-mapping
60+
i_line = find_i_line(lines, "variables:")
61+
lines[i_line + 1 : i_line + 1] = BASIC_CRS_STR.split("\n")
62+
63+
i_line = find_i_line(lines, "float node_data(")
64+
ref_lines = [
65+
# NOTE: space to match the surrounding indent
66+
f' node_data:grid_mapping = "crs: {includes}" ;',
67+
# NOTE: this is already present
68+
# f'node_data:coordinates = "{includes}" ;'
69+
]
70+
lines[i_line + 1 : i_line + 1] = ref_lines
71+
72+
cdl_str = "\n".join(lines)
73+
ncgen_from_cdl(cdl_str=cdl_str, cdl_path=None, nc_path=nc_path)
74+
75+
76+
@pytest.mark.parametrize("cs_axes", ["--", "x-", "-y", "xy"])
77+
def test_default_mesh_cs(tmp_path, cs_axes):
78+
"""Test coord-systems of mesh cube and coords, if location coords have a crs."""
79+
nc_path = tmp_path / "test_temp.nc"
80+
do_x = "x" in cs_axes
81+
do_y = "y" in cs_axes
82+
make_file(nc_path, node_x_crs=do_x, node_y_crs=do_y)
83+
with PARSE_UGRID_ON_LOAD.context():
84+
cube = iris.load_cube(nc_path, "node_data")
85+
meshco_x, meshco_y = [cube.coord(mesh_coords=True, axis=ax) for ax in ("x", "y")]
86+
# NOTE: at present, none of these load with a coordinate system,
87+
# because we don't support the extended grid-mapping syntax.
88+
# see: https://github.com/SciTools/iris/issues/3388
89+
assert meshco_x.coord_system is None
90+
assert meshco_y.coord_system is None
91+
92+
93+
def test_assigned_mesh_cs(tmp_path):
94+
# Check that when a coord system is manually assigned to a location coord,
95+
# the corresponding meshcoord reports the same cs.
96+
nc_path = tmp_path / "test_temp.nc"
97+
make_file(nc_path)
98+
with PARSE_UGRID_ON_LOAD.context():
99+
cube = iris.load_cube(nc_path, "node_data")
100+
nodeco_x = cube.mesh.coord(include_nodes=True, axis="x")
101+
meshco_x, meshco_y = [cube.coord(axis=ax) for ax in ("x", "y")]
102+
assert nodeco_x.coord_system is None
103+
assert meshco_x.coord_system is None
104+
assert meshco_y.coord_system is None
105+
assigned_cs = GeogCS(1.0)
106+
nodeco_x.coord_system = assigned_cs
107+
assert meshco_x.coord_system is assigned_cs
108+
assert meshco_y.coord_system is None
109+
# This also affects cube.coord_system(), even though it is an auxcoord,
110+
# since there are no dim-coords, or any other coord with a c-s.
111+
# TODO: this may be a mistake -- see https://github.com/SciTools/iris/issues/6051
112+
assert cube.coord_system() is assigned_cs
113+
114+
115+
def test_meshcoord_coordsys_copy(tmp_path):
116+
# Check that copying a meshcoord with a coord system works properly.
117+
nc_path = tmp_path / "test_temp.nc"
118+
make_file(nc_path)
119+
with PARSE_UGRID_ON_LOAD.context():
120+
cube = iris.load_cube(nc_path, "node_data")
121+
node_coord = cube.mesh.coord(include_nodes=True, axis="x")
122+
assigned_cs = GeogCS(1.0)
123+
node_coord.coord_system = assigned_cs
124+
mesh_coord = cube.coord(axis="x")
125+
assert mesh_coord.coord_system is assigned_cs
126+
meshco_copy = mesh_coord.copy()
127+
assert meshco_copy == mesh_coord
128+
# Note: still the same object, because it is derived from the same node_coord
129+
assert meshco_copy.coord_system is assigned_cs

0 commit comments

Comments
 (0)