Skip to content

Commit e962419

Browse files
trexfeatherstrexfeathers
authored
Merge edc82bb into d4489d1 
2 parents d4489d1 + edc82bb commit e962419

File tree

19 files changed

+1764
-328
lines changed

19 files changed

+1764
-328
lines changed

‎docs/src/further_topics/netcdf_io.rst‎

Lines changed: 157 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -235,3 +235,160 @@ Worked example:
235235
>>> my_coord.ignore_axis = True
236236
>>> print(guess_coord_axis(my_coord))
237237
None
238+
239+
Multiple Coordinate Systems and Ordered Axes
240+
--------------------------------------------
241+
242+
In a CF compliant NetCDF file, the coordinate variables associated with a
243+
data variable can specify a specific *coordinate system* that defines how
244+
the coordinate values relate to physical locations on the globe. For example,
245+
a coordinate might have values with units of metres that should be referenced
246+
against a *Transverse Mercator* projection with a specific origin. This
247+
information is not stored on the coordinate itself, but in a separate
248+
*grid mapping* variable. Furthermore, the grid mapping for a set of
249+
coordinates is associated with the data variable (not the coordinates
250+
variables) via the ``grid_mapping`` attribute.
251+
252+
For example, a temperature variable defined on a *rotated pole* grid might
253+
look like this in a NetCDF file (extract of relevant variables):
254+
255+
.. code-block:: text
256+
257+
float T(rlat,rlon) ;
258+
T:long_name = "temperature" ;
259+
T:units = "K" ;
260+
T:grid_mapping = "rotated_pole" ;
261+
262+
char rotated_pole ;
263+
rotated_pole:grid_mapping_name = "rotated_latitude_longitude" ;
264+
rotated_pole:grid_north_pole_latitude = 32.5 ;
265+
rotated_pole:grid_north_pole_longitude = 170. ;
266+
267+
float rlon(rlon) ;
268+
rlon:long_name = "longitude in rotated pole grid" ;
269+
rlon:units = "degrees" ;
270+
rlon:standard_name = "grid_longitude";
271+
272+
float rlat(rlat) ;
273+
rlat:long_name = "latitude in rotated pole grid" ;
274+
rlat:units = "degrees" ;
275+
rlat:standard_name = "grid_latitude";
276+
277+
278+
Note how the ``rotated pole`` grid mapping (coordinate system) is referenced
279+
from the data variable ``T:grid_mapping = "rotated_pole"`` and is implicitly
280+
associated with the dimension coordinate variables ``rlat`` and ``rlon``.
281+
282+
283+
Since version `1.8 of the CF Conventions
284+
<https://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.html#grid-mappings-and-projections>`_
285+
, there has been support for a more explicit version of the ``grid_mapping``
286+
attribute. This allows for **multiple coordinate systems** to be defined for
287+
a data variable and individual coordinates to be explicitly associated with
288+
a coordinate system. This is achieved by use of an **extended syntax** in the
289+
``grid_mapping`` variable of a data variable:
290+
291+
292+
.. code-block:: text
293+
294+
<grid_mapping_var>: <coord_var> [<coord_var>] [<grid_mapping_var>: <coord_var> ...]
295+
296+
where each ``grid_mapping_var`` identifies a grid mapping variable followed by
297+
the list of associated coordinate variables (``coord_var``). Note that with
298+
this syntax it is possible to specify multiple coordinate systems for a
299+
data variable.
300+
301+
For example, consider the following *air pressure* variable that is
302+
defined on an *OSGB Transverse Mercator grid*:
303+
304+
.. code-block:: text
305+
306+
float pres(y, x) ;
307+
pres:standard_name = "air_pressure" ;
308+
pres:units = "Pa" ;
309+
pres:coordinates = "lat lon" ;
310+
pres:grid_mapping = "crsOSGB: x y crsWGS84: lat lon" ;
311+
312+
double x(x) ;
313+
x:standard_name = "projection_x_coordinate" ;
314+
x:units = "m" ;
315+
316+
double y(y) ;
317+
y:standard_name = "projection_y_coordinate" ;
318+
y:units = "m" ;
319+
320+
double lat(y, x) ;
321+
lat:standard_name = "latitude" ;
322+
lat:units = "degrees_north" ;
323+
324+
double lon(y, x) ;
325+
lon:standard_name = "longitude" ;
326+
lon:units = "degrees_east" ;
327+
328+
int crsOSGB ;
329+
crsOSGB:grid_mapping_name = "transverse_mercator" ;
330+
crsOSGB:semi_major_axis = 6377563.396 ;
331+
crsOSGB:inverse_flattening = 299.3249646 ;
332+
<snip>
333+
334+
int crsWGS84 ;
335+
crsWGS84:grid_mapping_name = "latitude_longitude" ;
336+
crsWGS84:longitude_of_prime_meridian = 0. ;
337+
<snip>
338+
339+
340+
The dimension coordinates ``x`` and ``y`` are explicitly defined on
341+
an a *transverse mercator* grid via the ``crsOSGB`` variable.
342+
343+
However, with the extended grid syntax, it is also possible to define
344+
a second coordinate system on a standard **latitude_longitude** grid
345+
and associate it with the auxiliary ``lat`` and ``lon`` coordinates:
346+
347+
::
348+
349+
pres:grid_mapping = "crsOSGB: x y crsWGS84: lat lon" ;
350+
351+
352+
Note, the *order* of the axes in the extended grid mapping specification is
353+
significant, but only when used in conjunction with a
354+
`CRS Well Known Text (WKT)`_ representation of the coordinate system where it
355+
should be consistent with the ``AXES ORDER`` specified in the ``crs_wkt``
356+
attribute.
357+
358+
359+
Effect on loading
360+
^^^^^^^^^^^^^^^^^
361+
362+
When Iris loads a NetCDF file that uses the extended grid mapping syntax
363+
it will generate an :class:`iris.coord_systems.CoordSystem` for each
364+
coordinate system listed and attempt to attach it to the associated
365+
:class:`iris.coords.Coord` instances on the cube. Currently, Iris considers
366+
the ``crs_wkt`` supplementary and builds coordinate systems exclusively
367+
from the ``grid_mapping`` attribute.
368+
369+
The :attr:`iris.cube.Cube.extended_grid_mapping` property will be set to
370+
``True`` for cubes loaded from NetCDF data variables utilising the extended
371+
``grid_mapping`` syntax.
372+
373+
Effect on saving
374+
^^^^^^^^^^^^^^^^
375+
376+
To maintain existing behaviour, saving an :class:`iris.cube.Cube` to
377+
a netCDF file will default to the "simple" grid mapping syntax, unless
378+
the cube was loaded from a file using the extended grid mapping syntax.
379+
If the cube contains multiple coordinate systems, only the coordinate
380+
system of the dimension coordinate(s) will be specified.
381+
382+
To enable saving of multiple coordinate systems with ordered axes,
383+
set the :attr:`iris.cube.Cube.extended_grid_mapping` to ``True``.
384+
This will generate a ``grid_mapping`` attribute using the extended syntax
385+
to specify all coordinate systems on the cube. The axes ordering of the
386+
associated coordinate variables will be consistent with that of the
387+
generated ``crs_wkt`` attribute.
388+
389+
Note, the ``crs_wkt`` attribute will only be generated when the
390+
extended grid mapping is also written, i.e. when
391+
``Cube.extended_grid_mapping=True``.
392+
393+
394+
.. _CRS Well Known Text (WKT): https://cfconventions.org/Data/cf-conventions/cf-conventions-1.12/cf-conventions.html#use-of-the-crs-well-known-text-format

‎docs/src/whatsnew/latest.rst‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,31 @@ This document explains the changes made to Iris for this release
4040
now preserved when it was not previously. See also: :ref:`load-problems`.
4141
(:pull:`6465`, :pull:`6529`)
4242

43+
#. `@wjbenfold`_ and `@trexfeathers`_ added ``crs_wkt`` to the attributes when
44+
saving a :class:`~iris.coord_systems.CoordSystem` to a NetCDF file. Note that
45+
``crs_wkt`` is considered *supplementary* by the CF conventions, with
46+
``grid_mapping`` being the primary source of information, and ``crs_wkt`` not
47+
expected to contain conflicting information. Because of this, Iris generates
48+
:class:`~iris.coord_systems.CoordSystem` exclusively from ``grid_mapping``
49+
when loading, and writes a fresh ``crs_wkt`` whenever a
50+
:class:`~iris.coord_systems.CoordSystem` is saved. If your use case goes
51+
beyond the CF conventions, you can modify the save and load process for your
52+
needs by using the `Ncdata`_ package.
53+
See `CRS WKT in the CF Conventions`_ for more. (:issue:`3796`, :pull:`6519`)
54+
55+
#. `@ukmo-ccbunney`_ and `@trexfeathers`_ added support for
56+
**multiple coordinate systems** and **ordered coordinates** when loading
57+
and saving NetCDF files.
58+
This allows for coordinates to be explicitly associated with a coordinate
59+
system via an extended syntax in the ``grid_mapping`` attribute of a NetCDF
60+
data variable. This extended syntax also supports specification of multiple
61+
coordinate systems per data variable. Setting the property
62+
``cube.extended_grid_mapping = True`` will enable extended grid mapping
63+
syntax when saving a NetCDF file and also generate an associated **well known
64+
text** attribute (``crs_wkt``; as described in :issue:`3796`).
65+
See `CRS Grid Mappings and Projections`_ for more information.
66+
(:issue:`3388`:, :pull:`6536`:)
67+
4368
#. `@ESadek-MO`_ made MeshCoords immutable. :class:`iris.MeshCoord`s are now updated automatically when
4469
changing the attached mesh. All changes to the :class:`iris.MeshCoord` should instead be done to
4570
the relevant :class:`iris.Coord` located on the attached :class:`iris.MeshXY`. This change also affects
@@ -166,4 +191,7 @@ This document explains the changes made to Iris for this release
166191
.. comment
167192
Whatsnew resources in alphabetical order:
168193
194+
.. _CRS WKT in the CF Conventions: https://cfconventions.org/Data/cf-conventions/cf-conventions-1.12/cf-conventions.html#use-of-the-crs-well-known-text-format
195+
.. _CRS Grid Mappings and Projections: https://cfconventions.org/Data/cf-conventions/cf-conventions-1.12/cf-conventions.html#grid-mappings-and-projections
196+
.. _Ncdata: https://github.com/pp-mo/ncdata
169197
.. _Trusted Publishing: https://docs.pypi.org/trusted-publishers/

‎lib/iris/cube.py‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2461,6 +2461,32 @@ def coord_system(
24612461

24622462
return result
24632463

2464+
def coord_systems(self) -> list[iris.coord_systems.CoordSystem]:
2465+
"""Return a list of all coordinate systems used in cube coordinates."""
2466+
# Gather list of our unique CoordSystems on cube:
2467+
coord_systems = ClassDict(iris.coord_systems.CoordSystem)
2468+
for coord in self.coords():
2469+
if coord.coord_system:
2470+
coord_systems.add(coord.coord_system, replace=True)
2471+
2472+
return list(coord_systems.values())
2473+
2474+
@property
2475+
def extended_grid_mapping(self) -> bool:
2476+
"""Return True if a cube will use extended grid mapping syntax to write axes order in grid_mapping.
2477+
2478+
Only relevant when saving a cube to NetCDF file format.
2479+
2480+
For more details see "Grid Mappings and Projections" in the CF Conventions document:
2481+
https://cfconventions.org/cf-conventions/conformance.html
2482+
"""
2483+
return self.attributes.get("iris_extended_grid_mapping", False)
2484+
2485+
@extended_grid_mapping.setter
2486+
def extended_grid_mapping(self, ordered: bool) -> None:
2487+
"""Set to True to enable extended grid mapping syntax."""
2488+
self.attributes["iris_extended_grid_mapping"] = ordered
2489+
24642490
def _any_meshcoord(self) -> MeshCoord | None:
24652491
"""Return a MeshCoord if there are any, else None."""
24662492
mesh_coords = self.coords(mesh_coords=True)

‎lib/iris/exceptions.py‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -172,3 +172,9 @@ def __str__(self):
172172
"operations are not currently supported."
173173
)
174174
return msg.format(super().__str__())
175+
176+
177+
class CFParseError(IrisError):
178+
"""Raised when a string associated with a CF defined syntax could not be parsed."""
179+
180+
pass

‎lib/iris/fileformats/_nc_load_rules/actions.py‎

Lines changed: 65 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -203,7 +203,10 @@ def action_provides_grid_mapping(engine, gridmapping_fact):
203203

204204
def build_outer(engine_, cf_var_):
205205
coordinate_system = builder(engine_, cf_var_)
206-
engine_.cube_parts["coordinate_system"] = coordinate_system
206+
# We can now handle more than one coordinate_system, so store as dictionary:
207+
engine_.cube_parts["coordinate_systems"][cf_var_.cf_name] = (
208+
coordinate_system
209+
)
207210

208211
# Part 1 - only building - adding takes place downstream in
209212
# helpers.build_and_add_dimension/auxiliary_coordinate().
@@ -218,11 +221,8 @@ def build_outer(engine_, cf_var_):
218221
),
219222
)
220223

221-
# Check there is not an existing one.
222-
# ATM this is guaranteed by the caller, "run_actions".
223-
assert engine.fact_list("grid-type") == []
224-
225-
engine.add_fact("grid-type", (grid_mapping_type,))
224+
# Store grid-mapping name along with grid-type to match them later on
225+
engine.add_fact("grid-type", (var_name, grid_mapping_type))
226226

227227
else:
228228
message = "Coordinate system not created. Debug info:\n"
@@ -343,7 +343,35 @@ def action_build_dimension_coordinate(engine, providescoord_fact):
343343
# Non-conforming lon/lat/projection coords will be classed as
344344
# dim-coords by cf.py, but 'action_provides_coordinate' will give them
345345
# a coord-type of 'miscellaneous' : hence, they have no coord-system.
346-
coord_system = engine.cube_parts.get("coordinate_system")
346+
#
347+
# At this point, we need to match any "coordinate_system" entries in
348+
# the engine to the coord we are building. There are a couple of cases here:
349+
# 1. Simple `grid_mapping = crs` is used, in which case
350+
# we should just apply that mapping to all dim coords.
351+
# 2. Extended `grid_mapping = crs: coord1 coord2 crs: coord3 coord4`
352+
# is used in which case we need to match the crs to the coord here.
353+
354+
# We can have multiple coordinate_system, so now stored as a list (note plural key)
355+
coord_systems = engine.cube_parts.get("coordinate_systems")
356+
coord_system = None
357+
358+
if len(coord_systems):
359+
# Find which coord system applies to this coordinate.
360+
cs_mappings = engine.cube_parts.get("coordinate_system_mappings")
361+
if cs_mappings and coord_systems:
362+
if len(coord_systems) == 1 and None in cs_mappings:
363+
# Simple grid mapping (a single coord_system with no explicit coords)
364+
# Applies to spatial DimCoord(s) only. In this case only one
365+
# coordinate_system will have been built, so just use it.
366+
(coord_system,) = coord_systems.values()
367+
(cs_name,) = cs_mappings.values()
368+
else:
369+
# Extended grid mapping, e.g.
370+
# `grid_mapping = "crs: coord1 coord2 crs: coord3 coord4"`
371+
# We need to search for coord system that references our coordinate.
372+
if cs_name := cs_mappings.get(cf_var.cf_name):
373+
coord_system = coord_systems.get(cs_name, None)
374+
347375
# Translate the specific grid-mapping type to a grid-class
348376
if coord_system is None:
349377
succeed = True
@@ -352,8 +380,13 @@ def action_build_dimension_coordinate(engine, providescoord_fact):
352380
# Get a grid-class from the grid-type
353381
# i.e. one of latlon/rotated/projected, as for coord_grid_class.
354382
gridtypes_factlist = engine.fact_list("grid-type")
355-
(gridtypes_fact,) = gridtypes_factlist # only 1 fact
356-
(cs_gridtype,) = gridtypes_fact # fact contains 1 term
383+
384+
# potentially multiple grid-type facts; find one for CRS varname
385+
cs_gridtype = None
386+
for fact_cs_name, fact_cs_type in gridtypes_factlist:
387+
if fact_cs_name == cs_name:
388+
cs_gridtype = fact_cs_type
389+
357390
if cs_gridtype == "latitude_longitude":
358391
cs_gridclass = "latlon"
359392
elif cs_gridtype == "rotated_latitude_longitude":
@@ -446,6 +479,7 @@ def action_build_auxiliary_coordinate(engine, auxcoord_fact):
446479
"""Convert a CFAuxiliaryCoordinateVariable into a cube aux-coord."""
447480
(var_name,) = auxcoord_fact
448481
rule_name = "fc_build_auxiliary_coordinate"
482+
cf_var = engine.cf_var.cf_group[var_name]
449483

450484
# Identify any known coord "type" : latitude/longitude/time/time_period
451485
# If latitude/longitude, this sets the standard_name of the built AuxCoord
@@ -473,8 +507,27 @@ def action_build_auxiliary_coordinate(engine, auxcoord_fact):
473507
if coord_type:
474508
rule_name += f"_{coord_type}"
475509

510+
# Check if we have a coord_system specified for this coordinate.
511+
# (Only possible via extended grid_mapping attribute)
512+
coord_systems = engine.cube_parts.get("coordinate_systems")
513+
coord_system = None
514+
515+
cs_mappings = engine.cube_parts.get("coordinate_system_mappings", None)
516+
if cs_mappings and coord_systems:
517+
if len(coord_systems) == 1 and None in cs_mappings:
518+
# Simple grid_mapping - doesn't apply to AuxCoords (we need an explicit mapping)
519+
pass
520+
else:
521+
# Extended grid mapping, e.g.
522+
# `grid_mapping = "crs: coord1 coord2 crs: coord3 coord4"`
523+
# We need to search for coord system that references our coordinate.
524+
if cs_name := cs_mappings.get(cf_var.cf_name):
525+
coord_system = coord_systems.get(cs_name, None)
526+
476527
cf_var = engine.cf_var.cf_group.auxiliary_coordinates[var_name]
477-
hh.build_and_add_auxiliary_coordinate(engine, cf_var, coord_name=coord_name)
528+
hh.build_and_add_auxiliary_coordinate(
529+
engine, cf_var, coord_name=coord_name, coord_system=coord_system
530+
)
478531

479532
return rule_name
480533

@@ -615,10 +668,9 @@ def run_actions(engine):
615668
# default (all cubes) action, always runs
616669
action_default(engine) # This should run the default rules.
617670

618-
# deal with grid-mappings
671+
# deal with grid-mappings; potentially multiple mappings if extended grid_mapping used.
619672
grid_mapping_facts = engine.fact_list("grid_mapping")
620-
# For now, there should be at most *one* of these.
621-
assert len(grid_mapping_facts) in (0, 1)
673+
622674
for grid_mapping_fact in grid_mapping_facts:
623675
action_provides_grid_mapping(engine, grid_mapping_fact)
624676

0 commit comments

Comments
 (0)