Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Mobt 496 enforce forecast consistency #1900

Merged
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 0 additions & 68 deletions improver/calibration/reliability_calibration.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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
110 changes: 110 additions & 0 deletions improver/cli/enforce_consistent_forecasts.py
Original file line number Diff line number Diff line change
@@ -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
brhooper marked this conversation as resolved.
Show resolved Hide resolved
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)
160 changes: 160 additions & 0 deletions improver/utilities/enforce_consistency.py
Original file line number Diff line number Diff line change
@@ -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
brhooper marked this conversation as resolved.
Show resolved Hide resolved
# 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
12 changes: 9 additions & 3 deletions improver_tests/acceptance/SHA256SUMS
Original file line number Diff line number Diff line change
Expand Up @@ -170,9 +170,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
Expand Down
Loading