Support providing a static additional predictor when applying EMOS coefficients

improver/calibration/__init__.py

@@ -121,22 +121,24 @@ def split_forecasts_and_truth(
 def split_forecasts_and_coeffs(
     cubes: CubeList, land_sea_mask_name: Optional[str] = None,
 ):
-    """Split the input forecast, coefficients, land sea-mask and
-    probability template, if provided. The coefficients cubes and
-    land-sea mask are identified based on their name. The current
-    forecast and probability template are then split.
+    """Split the input forecast, coefficients, static additional predictors,
+    land sea-mask and probability template, if provided. The coefficients
+    cubes and land-sea mask are identified based on their name. The
+    static additional predictors are identified as not have a time
+    coordinate. The current forecast and probability template are then split.
 
     Args:
         cubes:
             A list of input cubes which will be split into relevant groups.
-            This includes the forecast, coefficients, land-sea mask and
-            probability template.
+            This includes the forecast, coefficients, static additional
+            predictors, land-sea mask and probability template.
         land_sea_mask_name:
             Name of the land-sea mask cube to help identification.
 
     Returns:
         - A cube containing the current forecast.
         - If found, a cubelist containing the coefficients else None.
+        - If found, a cubelist containing the static additional predictor else None.
         - If found, a land-sea mask will be returned, else None.
         - If found, a probability template will be returned, else None.
 
@@ -147,13 +149,16 @@ def split_forecasts_and_coeffs(
     coefficients = CubeList()
     land_sea_mask = None
     grouped_cubes: Dict[str, List[Cube]] = {}
+    static_additional_predictors = CubeList()
 
     for cubelist in cubes:
         for cube in cubelist:
             if "emos_coefficient" in cube.name():
                 coefficients.append(cube)
             elif land_sea_mask_name and cube.name() == land_sea_mask_name:
                 land_sea_mask = cube
+            elif "time" not in [c.name() for c in cube.coords()]:
+                static_additional_predictors.append(cube)
             else:
                 if "probability" in cube.name() and any(
                     "probability" in k for k in grouped_cubes
@@ -190,9 +195,13 @@ def split_forecasts_and_coeffs(
                 (current_forecast,) = grouped_cubes[key]
 
     coefficients = coefficients if coefficients else None
+    static_additional_predictors = (
+        static_additional_predictors if static_additional_predictors else None
+    )
     return (
         current_forecast,
         coefficients,
+        static_additional_predictors,
         land_sea_mask,
         prob_template,
     )

improver/calibration/ensemble_calibration.py

@@ -1721,6 +1721,45 @@ def __init__(self, percentiles: Optional[Sequence] = None):
         """
         self.percentiles = [np.float32(p) for p in percentiles] if percentiles else None
 
+    def _check_additional_field_sites(self, forecast, additional_fields):
+        """Check that the forecast and additional fields have matching sites.
+
+        Args:
+            forecast:
+                Uncalibrated forecast as probabilities, percentiles or
+                realizations
+            additional_fields:
+                Additional fields to be used as forecast predictors.
+
+        Raises:
+            ValueError: If the sites mismatch between the forecast and
+                additional fields.
+        """
+        if additional_fields:
+            if any([c.name() == "wmo_id" for c in forecast.coords()]):
+                sites = [
+                    np.array_equal(
+                        p.coord("wmo_id").points, forecast.coord("wmo_id").points
+                    )
+                    for p in additional_fields
+                ]
+                if not np.all(sites):
+                    mismatching_sites = []
+                    for ap in additional_fields:
+                        mismatching_sites.extend(
+                            list(
+                                set(ap.coord("wmo_id").points).symmetric_difference(
+                                    set(forecast.coord("wmo_id").points)
+                                )
+                            )
+                        )
+                    msg = (
+                        "The forecast and additional predictors have "
+                        f"mismatching sites. The mismatching sites are: "
+                        f"{list(set(mismatching_sites))}"
+                    )
+                    raise ValueError(msg)
+
     @staticmethod
     def _get_attribute(
         coefficients: CubeList, attribute_name: str, optional: bool = False
@@ -1958,6 +1997,7 @@ def process(
                 f"was {self.output_forecast_type}."
             )
             raise ValueError(msg)
+        self._check_additional_field_sites(forecast, additional_fields)
 
         forecast_as_realizations = forecast.copy()
         if self.input_forecast_type != "realizations":

improver/cli/apply_emos_coefficients.py

@@ -64,6 +64,9 @@ def process(
             could be either realizations, probabilities or percentiles.
             - A cubelist containing the coefficients used for calibration or None.
             If none then then input is returned unchanged.
+            - Optionally, cubes representing static additional predictors.
+            These static additional predictors is expected not to have a
+            time coordinate.
             - Optionally, a cube containing the land-sea mask on the same domain
             as the forecast that is to be calibrated. Land points are
             specified by ones and sea points are specified by zeros.
@@ -130,9 +133,13 @@ def process(
     from improver.calibration import split_forecasts_and_coeffs
     from improver.calibration.ensemble_calibration import ApplyEMOS
 
-    (forecast, coefficients, land_sea_mask, prob_template) = split_forecasts_and_coeffs(
-        cubes, land_sea_mask_name
-    )
+    (
+        forecast,
+        coefficients,
+        additional_predictors,
+        land_sea_mask,
+        prob_template,
+    ) = split_forecasts_and_coeffs(cubes, land_sea_mask_name)
 
     if coefficients is None:
         msg = (
@@ -146,6 +153,7 @@ def process(
     result = calibration_plugin(
         forecast,
         coefficients,
+        additional_fields=additional_predictors,
         land_sea_mask=land_sea_mask,
         prob_template=prob_template,
         realizations_count=realizations_count,

improver/ensemble_copula_coupling/ensemble_copula_coupling.py

@@ -988,8 +988,9 @@ def __repr__(self) -> str:
     def _check_template_cube(self, cube: Cube) -> None:
         """
         The template cube is expected to contain a leading threshold dimension
-        followed by spatial (y/x) dimensions. This check raises an error if
-        this is not the case. If the cube contains the expected dimensions,
+        followed by spatial (y/x) dimensions for a gridded cube. For a spot
+        template cube, the spatial dimensions are not expected to be dimension
+        coordinates. If the cube contains the expected dimensions,
         a threshold leading order is enforced.
 
         Args:
@@ -1000,7 +1001,8 @@ def _check_template_cube(self, cube: Cube) -> None:
         Raises:
             ValueError: If cube is not of the expected dimensions.
         """
-        check_for_x_and_y_axes(cube, require_dim_coords=True)
+        require_dim_coords = False if cube.coords("wmo_id") else True
+        check_for_x_and_y_axes(cube, require_dim_coords=require_dim_coords)
         dim_coords = get_dim_coord_names(cube)
         msg = (
             "{} expects a cube with only a leading threshold dimension, "

improver/utilities/cube_checker.py

@@ -50,7 +50,7 @@ def check_for_x_and_y_axes(cube: Cube, require_dim_coords: bool = False) -> None
 
     Raises:
         ValueError : Raise an error if non-uniform increments exist between
-                      grid points.
+            grid points.
     """
     for axis in ["x", "y"]:
         if require_dim_coords:
@@ -61,9 +61,7 @@ def check_for_x_and_y_axes(cube: Cube, require_dim_coords: bool = False) -> None
         if coord:
             pass
         else:
-            msg = "The cube does not contain the expected {}" "coordinates.".format(
-                axis
-            )
+            msg = "The cube does not contain the expected {} coordinates.".format(axis)
             raise ValueError(msg)
 
 

improver_tests/acceptance/SHA256SUMS

@@ -17,12 +17,18 @@ bd56f09033f7a37323e04be7b509e8e363eb2474c04fb8aa4fb5e3769acb81f7  ./apply-emos-c
 12110777cf3db3c8d76daa06f78760f8e314038262f359e659d16fa5d49b98bf  ./apply-emos-coefficients/realizations/realizations_coefficients.nc
 7c5cd2382cad998e22ef0d4210bf65d131587514de2b5cf4ad1fe3c30fede188  ./apply-emos-coefficients/realizations/realizations_kgo.nc
 615ddb59a838024bb055b9cbbb44b7c18ce28d5734952777774980760e0f36fc  ./apply-emos-coefficients/rebadged_percentiles/input.nc
-0bc91af1e7003d696aa440a07d79e653bca566733cdf9bf525ca92cff821548f  ./apply-emos-coefficients/sites/input.nc
+9b48f3eabeed93a90654a5e6a2f06fcb7297cdba40f650551d99a5e310cf27dd  ./apply-emos-coefficients/sites/additional_predictor/altitude.nc
+84fc5fbdc782ecff9262f26232a971b327e3adbf3f0b122d8c3d4b147fec172b  ./apply-emos-coefficients/sites/additional_predictor/coefficients.nc
+e610e4329ef038a5f91e190a3312c24b3e742000a35068268a8c17b9bc83944e  ./apply-emos-coefficients/sites/additional_predictor/percentile_kgo.nc
+60edea7a92af49d6afd10c687fd3f780caa8cd6f8a3bc29d33d5b3e618c949d8  ./apply-emos-coefficients/sites/additional_predictor/probability_kgo.nc
+a3bfead4d08e7ec2ac66067cd16d041b1aabe32eb242f3e072c67efd1d10d8f2  ./apply-emos-coefficients/sites/additional_predictor/probability_template.nc
 a0112b5a48ba99a1a4a345d43b0c453caaf25181504fb9a13786b5722f84cc10  ./apply-emos-coefficients/sites/offset/coefficients.nc
 22d950494dc9d76d642b0ecc75414e95bfb78fa693e5701fe7cd2fc462a0c6a7  ./apply-emos-coefficients/sites/offset/kgo.nc
 2722b1d08a87cebf1f36ac914d5badb9af38859b2150da83311dc202ce9be4dd  ./apply-emos-coefficients/sites/offset/offset_input.nc
+3c496867a4a77a4627426b79d8d963a83590b5440984cc2d7a31085cc0efbe10  ./apply-emos-coefficients/sites/percentile_input.nc
 63291b3435066f0b0fc56bc78b9774e0e7457d5a6b20086a67ecc6eb826d618d  ./apply-emos-coefficients/sites/point_by_point/coefficients.nc
 833a703e14a02999a7d70e9a37ce199ab36a09b13b96f87295a0fb91739658ce  ./apply-emos-coefficients/sites/point_by_point/kgo.nc
+0bc91af1e7003d696aa440a07d79e653bca566733cdf9bf525ca92cff821548f  ./apply-emos-coefficients/sites/realization_input.nc
 965a1f0565f9192d95eb01d0a61dc7bede6902f5e1a0481d92611e677841a139  ./apply-emos-coefficients/truncated_normal/input.nc
 feb2d93cf96b05889571bcb348dfbb4d60cfc1f3d9b7343dcf7d85cf34339746  ./apply-emos-coefficients/truncated_normal/kgo.nc
 d2cab1d3d8aa588be08a3d7d65e95e859fed37daa767e8d4a2bdaae25702b9a8  ./apply-emos-coefficients/truncated_normal/truncated_normal_coefficients.nc

improver_tests/acceptance/test_apply_emos_coefficients.py

@@ -84,7 +84,7 @@ def test_normal_point_by_point_sites(tmp_path):
     independently at each site (initial guess and minimisation)."""
     kgo_dir = acc.kgo_root() / "apply-emos-coefficients/sites/point_by_point"
     kgo_path = kgo_dir / "kgo.nc"
-    input_path = kgo_dir / ".." / "input.nc"
+    input_path = kgo_dir / ".." / "realization_input.nc"
     emos_est_path = kgo_dir / "coefficients.nc"
     output_path = tmp_path / "output.nc"
     args = [
@@ -352,6 +352,57 @@ def test_percentiles_in_probabilities_out(tmp_path):
     acc.compare(output_path, kgo_path, atol=LOOSE_TOLERANCE, rtol=LOOSE_TOLERANCE)
 
 
+def test_percentile_sites_additional_predictor(tmp_path):
+    """Test using percentile site forecasts with a static additional
+    predictor."""
+    kgo_dir = acc.kgo_root() / "apply-emos-coefficients/sites/additional_predictor"
+    kgo_path = kgo_dir / "percentile_kgo.nc"
+    input_path = kgo_dir / ".." / "percentile_input.nc"
+    emos_est_path = kgo_dir / "coefficients.nc"
+    additional_predictor_path = kgo_dir / "altitude.nc"
+    output_path = tmp_path / "output.nc"
+    args = [
+        input_path,
+        emos_est_path,
+        additional_predictor_path,
+        "--realizations-count",
+        "19",
+        "--random-seed",
+        "0",
+        "--output",
+        output_path,
+    ]
+    run_cli(args)
+    acc.compare(output_path, kgo_path, atol=LOOSE_TOLERANCE)
+
+
+def test_perc_in_prob_out_sites_additional_predictor(tmp_path):
+    """Test using percentile site forecasts with a static additional
+    predictor and a probability template to generate a probability
+    site forecast."""
+    kgo_dir = acc.kgo_root() / "apply-emos-coefficients/sites/additional_predictor"
+    kgo_path = kgo_dir / "probability_kgo.nc"
+    input_path = kgo_dir / ".." / "percentile_input.nc"
+    emos_est_path = kgo_dir / "coefficients.nc"
+    additional_predictor_path = kgo_dir / "altitude.nc"
+    prob_template = kgo_dir / "probability_template.nc"
+    output_path = tmp_path / "output.nc"
+    args = [
+        input_path,
+        emos_est_path,
+        additional_predictor_path,
+        prob_template,
+        "--realizations-count",
+        "19",
+        "--random-seed",
+        "0",
+        "--output",
+        output_path,
+    ]
+    run_cli(args)
+    acc.compare(output_path, kgo_path, atol=LOOSE_TOLERANCE)
+
+
 def test_no_coefficients(tmp_path):
     """Test no coefficients provided"""
     kgo_dir = acc.kgo_root() / "apply-emos-coefficients/normal"

improver_tests/calibration/ensemble_calibration/test_ApplyEMOS.py

@@ -30,6 +30,7 @@
 # POSSIBILITY OF SUCH DAMAGE.
 """Unit tests for the `ensemble_calibration.ApplyEMOS` class."""
 
+import datetime
 import unittest
 from typing import Sequence, Union
 
@@ -39,7 +40,10 @@
 from iris.tests import IrisTest
 
 from improver.calibration.ensemble_calibration import ApplyEMOS
+from improver.metadata.constants.attributes import MANDATORY_ATTRIBUTE_DEFAULTS
+from improver.spotdata.build_spotdata_cube import build_spotdata_cube
 from improver.synthetic_data.set_up_test_cubes import (
+    construct_scalar_time_coords,
     set_up_percentile_cube,
     set_up_probability_cube,
     set_up_variable_cube,
@@ -186,6 +190,55 @@ def setUp(self):
             land_sea_data, name="land_binary_mask", units="1"
         )
 
+        # Generate site forecast and additional predictor cubes.
+        data = np.tile([1.6, 1.3, 1.4, 1.1], (4, 1))
+        altitude = np.array([10, 20, 30, 40])
+        latitude = np.linspace(58.0, 59.5, 4)
+        longitude = np.linspace(-0.25, 0.5, 4)
+        wmo_id = ["03001", "03002", "03003", "03004"]
+        time_coords = construct_scalar_time_coords(
+            datetime.datetime(2017, 11, 5, 4, 0),
+            None,
+            datetime.datetime(2017, 11, 5, 0, 0),
+        )
+        time_coords = [t[0] for t in time_coords]
+        realization_coord = [
+            iris.coords.DimCoord(np.arange(1, 5), standard_name="realization")
+        ]
+        self.realizations_spot_cube = build_spotdata_cube(
+            data,
+            "air_temperature",
+            "degC",
+            altitude,
+            latitude,
+            longitude,
+            wmo_id,
+            scalar_coords=time_coords,
+            additional_dims=realization_coord,
+        )
+
+        self.realizations_spot_cube.attributes.update(MANDATORY_ATTRIBUTE_DEFAULTS)
+
+        self.spot_altitude_cube = self.realizations_spot_cube[0].copy(
+            self.realizations_spot_cube.coord("altitude").points
+        )
+        self.spot_altitude_cube.rename("altitude")
+        self.spot_altitude_cube.units = "m"
+        for coord in [
+            "altitude",
+            "forecast_period",
+            "forecast_reference_time",
+            "realization",
+            "time",
+        ]:
+            self.spot_altitude_cube.remove_coord(coord)
+
+        self.spot_coefficients = build_coefficients_cubelist(
+            self.realizations_spot_cube,
+            [0, [0.9, 0.1], 0, 1],
+            CubeList([self.realizations_spot_cube, self.spot_altitude_cube]),
+        )
+
     def test_null_percentiles(self):
         """Test effect of "neutral" emos coefficients in percentile space
         (this is small but non-zero due to limited sampling of the
@@ -267,6 +320,7 @@ def test_error_realizations_count(self):
             ApplyEMOS()(self.percentiles, self.coefficients)
 
     def test_additional_predictor(self):
+        """Test providing an additional predictor."""
         altitude = set_up_variable_cube(
             np.ones((3, 3), dtype=np.float32), name="surface_altitude", units="m"
         )
@@ -292,6 +346,30 @@ def test_additional_predictor(self):
         )
         self.assertArrayAlmostEqual(result.data, expected_data)
 
+    def test_realizations_additional_predictor_at_sites(self):
+        """Test providing an additional predictor for site forecasts."""
+        expected_data = np.tile([2.44, 3.17, 4.26, 4.99], (4, 1))
+        result = ApplyEMOS()(
+            self.realizations_spot_cube,
+            self.spot_coefficients,
+            additional_fields=CubeList([self.spot_altitude_cube]),
+            realizations_count=3,
+        )
+        self.assertArrayAlmostEqual(result.data, expected_data)
+
+    def test_additional_predictor_site_mismatch(self):
+        """Test for a mismatch in sites between the forecast and
+        the additional predictor."""
+        spot_altitude_cube = self.spot_altitude_cube[1:]
+        msg = "The forecast and additional predictors"
+        with self.assertRaisesRegex(ValueError, msg):
+            ApplyEMOS()(
+                self.realizations_spot_cube,
+                self.spot_coefficients,
+                additional_fields=CubeList([spot_altitude_cube]),
+                realizations_count=3,
+            )
+
     def test_land_sea_mask(self):
         """Test that coefficients can be effectively applied to "land" points
         only"""