Merge pull request #983 from MetOffice/aifs_comp

New callbacks to standardise some dim names

src/CSET/operators/_utils.py

@@ -66,6 +66,42 @@ def get_cube_yxcoordname(cube: iris.cube.Cube) -> tuple[str, str]:
     return (y_coords[0], x_coords[0])
 
 
+def is_spatialdim(cube: iris.cube.Cube) -> bool:
+    """Determine whether a cube is has two spatial dimension coordinates.
+
+    If cube has both spatial dims, it will contain two unique coordinates
+    that explain space (latitude and longitude). The coordinates have to
+    be iterable/contain usable dimension data, as cubes may contain these
+    coordinates as scalar dimensions after being collapsed.
+
+    Arguments
+    ---------
+    cube: iris.cube.Cube
+        An iris cube which will be checked to see if it contains coordinate
+        names that match a pre-defined list of acceptable coordinate names.
+
+    Returns
+    -------
+    bool
+        If true, then the cube has a spatial projection and thus can be plotted
+        as a map.
+    """
+    # Acceptable horizontal coordinate names.
+    X_COORD_NAMES = ["longitude", "grid_longitude", "projection_x_coordinate", "x"]
+    Y_COORD_NAMES = ["latitude", "grid_latitude", "projection_y_coordinate", "y"]
+
+    # Get a list of coordinate names for the cube
+    coord_names = [coord.name() for coord in cube.dim_coords]
+    x_coords = [coord for coord in coord_names if coord in X_COORD_NAMES]
+    y_coords = [coord for coord in coord_names if coord in Y_COORD_NAMES]
+
+    # If there is one coordinate for both x and y direction return True.
+    if len(x_coords) == 1 and len(y_coords) == 1:
+        return True
+    else:
+        return False
+
+
 def is_transect(cube: iris.cube.Cube) -> bool:
     """Determine whether a cube is a transect.
 

src/CSET/operators/read.py

@@ -198,6 +198,8 @@ def callback(cube: iris.cube.Cube, field, filename: str):
         _lfric_normalise_callback(cube, field, filename)
         _lfric_time_coord_fix_callback(cube, field, filename)
         _longitude_fix_callback(cube, field, filename)
+        _fix_spatialcoord_name_callback(cube)
+        _fix_pressurecoord_name_callback(cube)
 
     return callback
 
@@ -333,6 +335,52 @@ def _longitude_fix_callback(cube: iris.cube.Cube, field, filename):
     return cube
 
 
+def _fix_spatialcoord_name_callback(cube: iris.cube.Cube):
+    """Check latitude and longitude coordinates name.
+
+    This is necessary as some models define their grid as 'grid_latitude' and 'grid_longitude'
+    and this means that recipes will fail - particularly if the user is comparing multiple models
+    where the spatial coordinate names differ.
+    """
+    import CSET.operators._utils as utils
+
+    # Check if cube is spatial.
+    if not utils.is_spatialdim(cube):
+        # Don't modify non-spatial cubes.
+        return cube
+
+    # Get spatial coords.
+    y_name, x_name = utils.get_cube_yxcoordname(cube)
+
+    if y_name in ["latitude"] and cube.coord(y_name).units in [
+        "degrees",
+        "degrees_north",
+        "degrees_south",
+    ]:
+        cube.coord(y_name).rename("grid_latitude")
+    if x_name in ["longitude"] and cube.coord(x_name).units in [
+        "degrees",
+        "degrees_west",
+        "degrees_east",
+    ]:
+        cube.coord(x_name).rename("grid_longitude")
+
+
+def _fix_pressurecoord_name_callback(cube: iris.cube.Cube):
+    """Rename pressure_level coordinate to pressure if it exists.
+
+    This problem was raised because the AIFS model data from ECMWF
+    defines the pressure coordinate with the name 'pressure_level' rather
+    than compliant CF coordinate names
+    """
+    # We only want to modify instances where the coordinate system is actually
+    # latitude/longitude, and not touch the cube if the coordinate system is say
+    # meters.
+    for coord in cube.dim_coords:
+        if coord.name() == "pressure_level":
+            coord.rename("pressure")
+
+
 def _check_input_files(input_path: Path | str, filename_pattern: str) -> Iterable[Path]:
     """Get an iterable of files to load, and check that they all exist.
 

tests/operators/test_read.py

@@ -240,3 +240,36 @@ def test_lfric_time_coord_fix_callback_no_time():
     cube = iris.cube.Cube([0, 0, 0], aux_coords_and_dims=[(length_coord, 0)])
     read._lfric_time_coord_fix_callback(cube, None, None)
     assert len(cube.coords("time")) == 0
+
+
+def test_pressurecoordfix_callback():
+    """Check that pressure_level is renamed to pressure if it exists."""
+    cube = iris.load_cube("tests/test_data/transect_test_umpl.nc")
+    cube.coord("pressure").rename("pressure_level")
+    read._fix_pressurecoord_name_callback(cube)
+    assert (
+        repr(cube.coords())
+        == "[<DimCoord: time / (hours since 1970-01-01 00:00:00)  [...]  shape(2,)>, <DimCoord: pressure / (hPa)  [ 100., 150., ..., 950., 1000.]  shape(16,)>, <DimCoord: latitude / (degrees)  [-10.98, -10.94, ..., -10.82, -10.78]  shape(6,)>, <DimCoord: longitude / (degrees)  [19.02, 19.06, ..., 19.18, 19.22]  shape(6,)>, <DimCoord: forecast_reference_time / (hours since 1970-01-01 00:00:00)  [...]>, <AuxCoord: forecast_period / (hours)  [15., 18.]  shape(2,)>]"
+    )
+
+
+def test_spatialcoordrename_callback():
+    """Check that spatial coord gets renamed if it is not grid_latitude."""
+    # This cube contains 'latitude' and 'longitude'
+    cube = iris.load_cube("tests/test_data/transect_test_umpl.nc")
+    read._fix_spatialcoord_name_callback(cube)
+    assert (
+        repr(cube.coords())
+        == "[<DimCoord: time / (hours since 1970-01-01 00:00:00)  [...]  shape(2,)>, <DimCoord: pressure / (hPa)  [ 100., 150., ..., 950., 1000.]  shape(16,)>, <DimCoord: grid_latitude / (degrees)  [-10.98, -10.94, ..., -10.82, -10.78]  shape(6,)>, <DimCoord: grid_longitude / (degrees)  [19.02, 19.06, ..., 19.18, 19.22]  shape(6,)>, <DimCoord: forecast_reference_time / (hours since 1970-01-01 00:00:00)  [...]>, <AuxCoord: forecast_period / (hours)  [15., 18.]  shape(2,)>]"
+    )
+
+
+def test_spatialcoordnotexist_callback():
+    """Check that spatial coord returns cube if cube does not contain spatial coordinates."""
+    cube = iris.load_cube("tests/test_data/transect_test_umpl.nc")
+    cube = cube[:, :, 0, 0]  # Remove spatial dimcoords
+    read._fix_spatialcoord_name_callback(cube)
+    assert (
+        repr(cube.coords())
+        == "[<DimCoord: time / (hours since 1970-01-01 00:00:00)  [...]  shape(2,)>, <DimCoord: pressure / (hPa)  [ 100., 150., ..., 950., 1000.]  shape(16,)>, <DimCoord: forecast_reference_time / (hours since 1970-01-01 00:00:00)  [...]>, <DimCoord: latitude / (degrees)  [-10.98]>, <DimCoord: longitude / (degrees)  [19.02]>, <AuxCoord: forecast_period / (hours)  [15., 18.]  shape(2,)>]"
+    )

tests/operators/test_utils.py

@@ -14,30 +14,31 @@
 
 """Tests for common operator functionality across CSET."""
 
+import iris
 import pytest
 
 import CSET.operators._utils as operator_utils
 
 
 def test_missing_coord_get_cube_yxcoordname_x(regrid_rectilinear_cube):
     """Missing X coordinate raises error."""
-    regrid_rectilinear_cube.remove_coord("longitude")
+    regrid_rectilinear_cube.remove_coord("grid_longitude")
     with pytest.raises(ValueError):
         operator_utils.get_cube_yxcoordname(regrid_rectilinear_cube)
 
 
 def test_missing_coord_get_cube_yxcoordname_y(regrid_rectilinear_cube):
     """Missing Y coordinate raises error."""
-    regrid_rectilinear_cube.remove_coord("longitude")
+    regrid_rectilinear_cube.remove_coord("grid_longitude")
     with pytest.raises(ValueError):
         operator_utils.get_cube_yxcoordname(regrid_rectilinear_cube)
 
 
 def test_get_cube_yxcoordname(regrid_rectilinear_cube):
     """Check that function returns tuple containing horizontal dimension names."""
     assert (operator_utils.get_cube_yxcoordname(regrid_rectilinear_cube)) == (
-        "latitude",
-        "longitude",
+        "grid_latitude",
+        "grid_longitude",
     )
 
 
@@ -58,3 +59,16 @@ def test_is_transect_correctcoord(transect_source_cube):
     # Retain only time and latitude coordinate, so it passes the first spatial coord test.
     transect_source_cube_slice = transect_source_cube[:, :, :, 0]
     assert operator_utils.is_transect(transect_source_cube_slice)
+
+
+def test_is_spatialdim_false():
+    """Check that is spatial test returns false if cube does not contain spatial coordinates."""
+    cube = iris.load_cube("tests/test_data/transect_test_umpl.nc")
+    cube = cube[:, :, 0, 0]  # Remove spatial dimcoords
+    assert not operator_utils.is_spatialdim(cube)
+
+
+def test_is_spatialdim_true():
+    """Check that is spatial test returns true if cube contains spatial coordinates."""
+    cube = iris.load_cube("tests/test_data/transect_test_umpl.nc")
+    assert operator_utils.is_spatialdim(cube)