Generate calibrated forecasts from EMOS with an alternative percentil…

…e set (metoppv#1587)

* Add capability to pass a comma separated list of percentiles to the apply EMOS CLI. These will be the percentiles output within the calibrated forecast, if these percentiles are provided. Otherwise, the percentiles output will match the percentiles input.

* Minor updates.

improver/calibration/ensemble_calibration.py

@@ -38,7 +38,7 @@
 
 """
 import warnings
-from typing import Any, Callable, Dict, List, Optional, Tuple, Union
+from typing import Any, Callable, Dict, List, Optional, Sequence, Tuple, Union
 
 import iris
 import numpy as np
@@ -1712,6 +1712,15 @@ class ApplyEMOS(PostProcessingPlugin):
     Class to calibrate an input forecast given EMOS coefficients
     """
 
+    def __init__(self, percentiles: Optional[Sequence] = None):
+        """Initialise class.
+
+        Args:
+            percentiles:
+                The set of percentiles used to create the calibrated forecast.
+        """
+        self.percentiles = [np.float32(p) for p in percentiles] if percentiles else None
+
     @staticmethod
     def _get_attribute(
         coefficients: CubeList, attribute_name: str, optional: bool = False
@@ -1818,11 +1827,12 @@ def _convert_to_realizations(
 
         return forecast_as_realizations
 
-    def _calibrate_forecast(
+    def _format_forecast(
         self, forecast: Cube, randomise: bool, random_seed: int
     ) -> Cube:
         """
-        Generate calibrated probability, percentile or realization output
+        Generate calibrated probability, percentile or realization output in
+        the desired format.
 
         Args:
             forecast:
@@ -1860,7 +1870,9 @@ def _calibrate_forecast(
                     self.distribution["location"],
                     self.distribution["scale"],
                     forecast,
-                    percentiles=perc_coord.points,
+                    percentiles=self.percentiles
+                    if self.percentiles
+                    else perc_coord.points,
                 )
             else:
                 no_of_percentiles = len(forecast.coord("realization").points)
@@ -1952,7 +1964,7 @@ def process(
             ),
         }
 
-        result = self._calibrate_forecast(forecast, randomise, random_seed)
+        result = self._format_forecast(forecast, randomise, random_seed)
 
         if land_sea_mask:
             # fill in masked sea points with uncalibrated data

improver/cli/apply_emos_coefficients.py

@@ -56,6 +56,7 @@ def process(
     random_seed: int = None,
     ignore_ecc_bounds=False,
     predictor="mean",
+    percentiles: cli.comma_separated_list = None,
 ):
     """Applying coefficients for Ensemble Model Output Statistics.
 
@@ -108,6 +109,8 @@ def process(
             the location parameter when estimating the EMOS coefficients.
             Currently the ensemble mean ("mean") and the ensemble
             realizations ("realizations") are supported as the predictors.
+        percentiles (List[float]):
+            The set of percentiles used to create the calibrated forecast.
 
     Returns:
         iris.cube.Cube:
@@ -139,7 +142,7 @@ def process(
         msg = "The land_sea_mask cube does not have the name 'land_binary_mask'"
         raise ValueError(msg)
 
-    calibration_plugin = ApplyEMOS()
+    calibration_plugin = ApplyEMOS(percentiles=percentiles)
     result = calibration_plugin(
         cube,
         coefficients,

improver_tests/acceptance/SHA256SUMS

@@ -2,6 +2,7 @@
 a6a84d0142796e4b9ca7bd3f0ad78586ea77684f5df02732fdda2ab54233cbb6  ./aggregate-reliability-tables/basic/multiple_tables_kgo.nc
 c0c38c5b1ba16fd5f7310b6d2ee0ab9d8b6bcdb3b23c8475193eefae1eb14a17  ./aggregate-reliability-tables/basic/reliability_table.nc
 9eeda326cdc66d93e591c87e9f7ceb17cea329397169fb0518f50da3bf4dc61f  ./aggregate-reliability-tables/basic/reliability_table_2.nc
+52b02f4c6e4adc604d4961ce87b4ada5985d07eac1d3c703cbab985de223228f  ./apply-emos-coefficients/alternative_percentiles/kgo.nc
 ff3a00a16fa94697d6e529a6f98384e4602f70a7d6ce26669c3947c8ac2e6f7e  ./apply-emos-coefficients/land_sea/landmask.nc
 bee125371a9db4c80fb7f8d1b0644a2710903aad27023b5939331f7ec5f3b617  ./apply-emos-coefficients/land_sea/percentiles_kgo.nc
 7f5cabc9b5de4c13aebf833d5372159edf1deafd19bba8721318b889a09b407f  ./apply-emos-coefficients/land_sea/probabilities_kgo.nc

improver_tests/acceptance/test_apply_emos_coefficients.py

@@ -258,6 +258,28 @@ def test_percentiles_input_land_sea(tmp_path):
     acc.compare(output_path, kgo_path, atol=LOOSE_TOLERANCE)
 
 
+def test_alternative_percentiles(tmp_path):
+    """Test using percentiles as input with an alternative set of
+    percentiles specified."""
+    kgo_dir = acc.kgo_root() / "apply-emos-coefficients/alternative_percentiles"
+    kgo_path = kgo_dir / "kgo.nc"
+    input_path = kgo_dir / "../percentiles/input.nc"
+    emos_est_path = kgo_dir / "../normal/normal_coefficients.nc"
+    output_path = tmp_path / "output.nc"
+    args = [
+        input_path,
+        emos_est_path,
+        "--realizations-count",
+        "18",
+        "--percentiles",
+        "25,50,75",
+        "--output",
+        output_path,
+    ]
+    run_cli(args)
+    acc.compare(output_path, kgo_path, atol=LOOSE_TOLERANCE, rtol=LOOSE_TOLERANCE)
+
+
 def test_percentiles_error(tmp_path):
     """Test using percentiles as input"""
     kgo_dir = acc.kgo_root() / "apply-emos-coefficients/percentiles"

improver_tests/calibration/ensemble_calibration/test_ApplyEMOS.py

@@ -179,6 +179,7 @@ def setUp(self):
                 np.full((3, 3), 10.534898),
             ]
         )
+        self.alternative_percentiles = [25.0, 50.0, 75.0]
 
     def test_null_percentiles(self):
         """Test effect of "neutral" emos coefficients in percentile space
@@ -375,6 +376,27 @@ def test_null_percentiles_frt_fp_mismatch(self):
             np.mean(result.data), self.null_percentiles_expected_mean
         )
 
+    def test_alternative_percentiles(self):
+        """Test that the calibrated forecast is at a specified set of
+        percentiles."""
+        result = ApplyEMOS(percentiles=self.alternative_percentiles)(
+            self.percentiles, self.coefficients, realizations_count=3
+        )
+        self.assertArrayEqual(
+            result.coord("percentile").points, self.alternative_percentiles
+        )
+
+    def test_alternative_string_percentiles(self):
+        """Test that the calibrated forecast is at a specified set of
+        percentiles where the input percentiles are strings."""
+        str_percentiles = list(map(str, self.alternative_percentiles))
+        result = ApplyEMOS(percentiles=str_percentiles)(
+            self.percentiles, self.coefficients, realizations_count=3
+        )
+        self.assertArrayEqual(
+            result.coord("percentile").points, self.alternative_percentiles
+        )
+
     def test_invalid_attribute(self):
         """Test that an exception is raised if multiple different distribution
         attributes are provided within the coefficients cubelist."""