Mobt 496 enforce forecast consistency (metoppv#1900)

* Remove obsolete enforce_consistent_probabilities code. Update checksums

* Add plugin and CLI for general enforcement of forecast consistency. Add corresponding unit and acceptance tests.

* Fix error in acceptance test.

* Corrections to formatting.

* Remove obsolete module import.

* Update docstrings.

* Minor change following review.

* Minor changes following review

* Updated checksums following reduction in acceptance test data size

* Update docstrings to include applications to realization data. Add unit and acceptance tests for realization data. Update checksum for new acceptance test data.

* Add limitation that only the identity function can be applied to probability reference forecasts.

* Formatting corrections

* Minor changes following review. Update checksums for change to acceptance test data.

improver/calibration/reliability_calibration.py

@@ -49,7 +49,6 @@
     create_unified_frt_coord,
     filter_non_matching_cubes,
 )
-from improver.cube_combiner import Combine
 from improver.metadata.probabilistic import (
     find_threshold_coordinate,
     probability_is_above_or_below,
@@ -1442,70 +1441,3 @@ def process(
         calibrated_forecast.data = calibrated_forecast.data.astype("float32")
 
         return calibrated_forecast
-
-
-class EnforceConsistentProbabilities(PostProcessingPlugin):
-    """Reduces the probabilities of a forecast to make it consistent
-    with another provided reference forecast, such that at any grid square the
-    probability of the forecast is less than the reference forecast. As an example, to
-    reduce the probability of a low cloud forecast to be less than the probability of the
-    reference forecast of total cloud.
-    """
-
-    def __init__(self, diff_for_warning: float = None) -> None:
-        """
-        Initialise class for enforcing probabilities between two forecasts.
-
-        Args:
-            diff_for_warning:
-                A float between 0 and 1. If assigned, the plugin will raise a warning
-                if the forecast probabilities are decreased by more than this value.
-        """
-        self.diff_for_warning = diff_for_warning
-
-    def process(self, forecast_cube: Cube, ref_forecast: Cube) -> Cube:
-        """
-        Lowers the probabilities of the forecast_cube to make it consistent with
-        the reference forecast. This is done by lowering the probabilities of the forecast cube,
-        to be equal to the reference forecast if they are higher than the reference forecast.
-        If the probability is already less than or equal to the reference forecast then the
-        probability is not altered.
-
-        Args:
-            forecast_cube:
-                A forecast cube of probabilities
-            ref_forecast:
-                A cube of probabilities that is used as an upper limit for
-                the forecast_cube probabilities. It must have the same dimensions
-                as the forecast_cube.
-
-        Returns:
-            A forecast cube of probabilities that has identical metadata to forecast_cube but
-            with reduced probabilities to make it consistent with ref_forecast
-        """
-
-        # re-name ref_forecast's threshold coordinate to ensure the coordinate names match for
-        # the combine plugin
-        forecast_threshold = find_threshold_coordinate(forecast_cube)
-        ref_forecast_threshold = find_threshold_coordinate(ref_forecast)
-
-        ref_forecast_threshold.standard_name = forecast_threshold.standard_name
-        ref_forecast_threshold.long_name = forecast_threshold.long_name
-
-        diff = Combine(operation="-")([ref_forecast, forecast_cube])
-
-        # clip the positive values such that only negative values or zeroes remain.
-        # These differences will be added to the forecast_cube probabilities.
-        diff.data = np.clip(diff.data, None, 0)
-
-        new_forecast = Combine(operation="+")([forecast_cube, diff])
-
-        largest_change = abs(np.amin(diff.data))
-        if self.diff_for_warning is not None and largest_change > self.diff_for_warning:
-            warnings.warn(
-                f"Inconsistency between forecast {forecast_cube.name} and {ref_forecast.name}"
-                f"is greater than {self.diff_for_warning}. Maximum absolute difference reported"
-                f"was {largest_change}"
-            )
-
-        return new_forecast

improver/cli/enforce_consistent_forecasts.py

@@ -0,0 +1,110 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# -----------------------------------------------------------------------------
+# (C) British Crown copyright. The Met Office.
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+#
+# * Redistributions of source code must retain the above copyright notice, this
+#   list of conditions and the following disclaimer.
+#
+# * Redistributions in binary form must reproduce the above copyright notice,
+#   this list of conditions and the following disclaimer in the documentation
+#   and/or other materials provided with the distribution.
+#
+# * Neither the name of the copyright holder nor the names of its
+#   contributors may be used to endorse or promote products derived from
+#   this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+# POSSIBILITY OF SUCH DAMAGE.
+"""CLI to enforce consistency between two forecasts."""
+
+from improver import cli
+
+
+@cli.clizefy
+@cli.with_output
+def process(
+    *cubes: cli.inputcubelist,
+    ref_name: str = None,
+    additive_amount: float = None,
+    multiplicative_amount: float = None,
+    comparison_operator: str = ">=",
+    diff_for_warning: float = None,
+):
+    """
+    Module to enforce that the values in the forecast cube are not less than or not
+    greater than a linear function of the corresponding values in the
+    reference forecast.
+
+    Args:
+        cubes (iris.cube.CubeList or list of iris.cube.Cube):
+            containing:
+                forecast_cube (iris.cube.Cube):
+                    Cube of forecasts to be updated by using the reference forecast to
+                    create a bound on the value of the forecasts.
+                ref_forecast (iris.cube.Cube)
+                    Cube of forecasts used to create a bound for the values in
+                    forecast_cube. It must be the same shape as forecast_cube but have
+                    a different name.
+        ref_name (str): Name of ref_forecast cube
+        additive_amount (float): The amount to be added to the reference forecast (in
+            the units of the reference forecast) prior to enforcing consistency between
+            the forecast and reference forecast. If both an additive_amount and
+            multiplicative_amount are specified then addition occurs after
+            multiplication. This option cannot be used for probability forecasts, if it
+            is then an error will be raised.
+        multiplicative_amount (float): The amount to multiply the reference forecast by
+            prior to enforcing consistency between the forecast and reference
+            forecast. If both an additive_amount and multiplicative_amount are
+            specified then addition occurs after multiplication. This option cannot be
+            used for probability forecasts, if it is then an error will be raised.
+        comparison_operator (str): Determines whether the forecast is enforced to be not
+            less than or not greater than the reference forecast. Valid choices are
+            ">=", for not less than, and "<=" for not greater than.
+        diff_for_warning (float): If assigned, the plugin will raise a warning if any
+            absolute change in forecast value is greater than this value.
+
+    Returns:
+        iris.cube.Cube:
+            A forecast cube with identical metadata to forecast but the forecasts are
+            enforced to be not less than or not greater than a linear function of
+            reference_forecast.
+    """
+    from iris.cube import CubeList
+
+    from improver.utilities.enforce_consistency import EnforceConsistentForecasts
+    from improver.utilities.flatten import flatten
+
+    cubes = flatten(cubes)
+
+    if len(cubes) != 2:
+        raise ValueError(
+            f"Exactly two cubes should be provided but received {len(cubes)}"
+        )
+
+    ref_forecast = CubeList(cubes).extract_cube(ref_name)
+    cubes.remove(ref_forecast)
+
+    forecast_cube = cubes[0]
+
+    plugin = EnforceConsistentForecasts(
+        additive_amount=additive_amount,
+        multiplicative_amount=multiplicative_amount,
+        comparison_operator=comparison_operator,
+        diff_for_warning=diff_for_warning,
+    )
+
+    return plugin(forecast_cube, ref_forecast)

improver/utilities/enforce_consistency.py

@@ -0,0 +1,160 @@
+# -*- coding: utf-8 -*-
+# -----------------------------------------------------------------------------
+# (C) British Crown copyright. The Met Office.
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+#
+# * Redistributions of source code must retain the above copyright notice, this
+#   list of conditions and the following disclaimer.
+#
+# * Redistributions in binary form must reproduce the above copyright notice,
+#   this list of conditions and the following disclaimer in the documentation
+#   and/or other materials provided with the distribution.
+#
+# * Neither the name of the copyright holder nor the names of its
+#   contributors may be used to endorse or promote products derived from
+#   this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+# POSSIBILITY OF SUCH DAMAGE.
+"""Provides utilities for enforcing consistency between forecasts."""
+import warnings
+from typing import Optional
+
+import numpy as np
+from iris.cube import Cube
+
+from improver import PostProcessingPlugin
+
+
+class EnforceConsistentForecasts(PostProcessingPlugin):
+    """Enforce that the forecasts provided are no less than or no greater than some
+    linear function of a reference forecast. For example, wind speed forecasts may be
+    provided as the reference forecast for wind gust forecasts with a requirement that
+    wind gusts are not less than 110% of the corresponding wind speed forecast. Wind
+    speed forecasts of 10 m/s will mean that the resulting wind gust forecast will be
+    at least 11 m/s.
+    """
+
+    def __init__(
+        self,
+        additive_amount: float = None,
+        multiplicative_amount: float = None,
+        comparison_operator: str = ">=",
+        diff_for_warning: Optional[float] = None,
+    ) -> None:
+        """
+        Initialise class for enforcing consistency between a forecast and a linear
+        function of a reference forecast.
+
+        Args:
+            additive_amount: The amount to be added to the reference forecast prior to
+                enforcing consistency between the forecast and reference forecast. If
+                both an additive_amount and multiplicative_amount are specified then
+                addition occurs after multiplication. This option cannot be used for
+                probability forecasts, if it is then an error will be raised.
+            multiplicative_amount: The amount to multiply the reference forecast by
+                prior to enforcing consistency between the forecast and reference
+                forecast. If both an additive_amount and multiplicative_amount are
+                specified then addition occurs after multiplication. This option cannot
+                be used for probability forecasts, if it is then an error will be raised.
+            comparison_operator: Determines whether the forecast is enforced to be not
+                less than or not greater than the reference forecast. Valid choices are
+                ">=", for not less than, and "<=" for not greater than.
+            diff_for_warning: If assigned, the plugin will raise a warning if any
+                absolute change in forecast value is greater than this value.
+        """
+
+        self.additive_amount = additive_amount
+        self.multiplicative_amount = multiplicative_amount
+        self.comparison_operator = comparison_operator
+        self.diff_for_warning = diff_for_warning
+
+    def process(self, forecast: Cube, reference_forecast: Cube) -> Cube:
+        """
+        Function to enforce that the values in the forecast cube are not less than or
+        not greater than a linear function of the corresponding values in
+        reference_forecast.
+
+        Args:
+            forecast: A forecast cube
+            reference_forecast: A reference forecast cube used to determine the bound
+                of the forecast cube.
+
+        Returns:
+            A forecast cube with identical metadata to forecast but the forecasts are
+            enforced to be not less than or not greater than a linear function of
+            reference_forecast.
+        """
+        # check forecast and reference units match
+        if forecast.units != reference_forecast.units:
+            msg = (
+                f"The units in the forecast and reference cubes do not match. The"
+                f"units of forecast were {forecast.units}, the units of "
+                f"reference_forecast were {reference_forecast.units}."
+            )
+            raise ValueError(msg)
+
+        # linear transformation cannot be applied to probability forecasts
+        if self.additive_amount is not None or self.multiplicative_amount is not None:
+            if forecast.name().startswith("probability_of"):
+                msg = (
+                    "For probability data, the additive_amount and multiplicative_amount"
+                    "inputs cannot be used. The input additive_amount was "
+                    f"{self.additive_amount}, the input multiplicative amount was "
+                    f"{self.multiplicative_amount}."
+                )
+                raise ValueError(msg)
+
+        # if additive_amount or multiplicative_amount not specified, then use identity
+        # function
+        if self.additive_amount is None:
+            self.additive_amount = 0
+        if self.multiplicative_amount is None:
+            self.multiplicative_amount = 1
+
+        # calculate forecast_bound by applying specified linear transformation to
+        # reference_forecast
+        forecast_bound = reference_forecast.copy()
+        forecast_bound.data = self.multiplicative_amount * forecast_bound.data
+        forecast_bound.data = self.additive_amount + forecast_bound.data
+
+        # assign forecast_bound to be either an upper or lower bound depending on input
+        # comparison_operator
+        lower_bound = None
+        upper_bound = None
+        if self.comparison_operator == ">=":
+            lower_bound = forecast_bound.data
+        elif self.comparison_operator == "<=":
+            upper_bound = forecast_bound.data
+        else:
+            msg = (
+                f"Comparison_operator must be either '>=' or '<=', not "
+                f"{self.comparison_operator}."
+            )
+            raise ValueError(msg)
+
+        new_forecast = forecast.copy()
+        new_forecast.data = np.clip(new_forecast.data, lower_bound, upper_bound)
+
+        diff = new_forecast.data - forecast.data
+        max_abs_diff = np.max(np.abs(diff))
+        if self.diff_for_warning is not None and max_abs_diff > self.diff_for_warning:
+            warnings.warn(
+                f"Inconsistency between forecast {forecast.name} and "
+                f"{reference_forecast.name} is greater than {self.diff_for_warning}. "
+                f"Maximum absolute difference reported was {max_abs_diff}"
+            )
+
+        return new_forecast

improver_tests/acceptance/SHA256SUMS

@@ -176,9 +176,15 @@ b946c7687cb9ed02a12a934429a31306004ad45214cf4b451468b077018c0911  ./convection-r
 d3efbc6014743793fafed512a8017a7b75e4f0ffa7fd202cd4f1b4a1086b2583  ./create-grid-with-halo/basic/kgo.nc
 fee00437131d2367dc317e5b0fff44e65e03371b8f096bf1ac0d4cc7253693c9  ./create-grid-with-halo/basic/source_grid.nc
 bf7e42be7897606682c3ecdaeb27bf3d3b6ab13a9a88b46c88ae6e92801c6245  ./create-grid-with-halo/halo_size/kgo.nc
-5474c4c2309dcc0991355007f20869eb6d1d234d721bdf37542cddb0570d7217  ./enforce-consistent-probabilities/hybrid_total_cloud.nc
-7e435c2db5e792b71bb5170140fbe1e37cc03d96bd88f73ce7196f9469ba3e13  ./enforce-consistent-probabilities/kgo.nc
-764833f16f1ed5939ba82130fa73961fdecacceebcf2b4cdb2970cb670e2ee3e  ./enforce-consistent-probabilities/total_cloud.nc
+51f9ff2c8e6cad54d04d6323c654ce25b812bea3ba6b0d85df21c19731c580fc  ./enforce-consistent-forecasts/percentile_forecast.nc
+5f7fa625272af0a51f6a59dd20a6f75cb53aa46299630783c32a77242c607313  ./enforce-consistent-forecasts/percentile_kgo.nc
+e210bf956dd3574eda3a64fdc3371ab16f85048ca120a3823e90d751d2325c46  ./enforce-consistent-forecasts/percentile_reference.nc
+5474c4c2309dcc0991355007f20869eb6d1d234d721bdf37542cddb0570d7217  ./enforce-consistent-forecasts/probability_forecast.nc
+7e435c2db5e792b71bb5170140fbe1e37cc03d96bd88f73ce7196f9469ba3e13  ./enforce-consistent-forecasts/probability_kgo.nc
+764833f16f1ed5939ba82130fa73961fdecacceebcf2b4cdb2970cb670e2ee3e  ./enforce-consistent-forecasts/probability_reference.nc
+40a05397b22e614bce592c1644a8443a27a94e70914cb482117a2fcc11ea8408  ./enforce-consistent-forecasts/realization_forecast.nc
+09eb8d42b1df18831660bfebef46a09436e66946a48488d8be89157e8e3666ef  ./enforce-consistent-forecasts/realization_kgo.nc
+9604ec25ad2256d41f824c8be957ae20158d75a7a8d632af249b2778d20934d0  ./enforce-consistent-forecasts/realization_reference.nc
 9b48f3eabeed93a90654a5e6a2f06fcb7297cdba40f650551d99a5e310cf27dd  ./estimate-emos-coefficients-from-table/altitude.nc
 93e0c05333bc93ca252027ee5e2da12ee93fac5bbff4bab00b9c334ad83141e2  ./estimate-emos-coefficients-from-table/forecast_table/_common_metadata
 037e1cfc55c703a30cdbc1dee78da47020a598683247dc29a8ef88e65b53dfb1  ./estimate-emos-coefficients-from-table/forecast_table/_metadata