Merge remote-tracking branch 'upstream/master' into improver1538_tabular_ingestion_functions

* upstream/master:
  Generate calibrated forecasts from EMOS with an alternative percentile set (metoppv#1587)
  MOBT-77: Weather symbols to represent an extended period (metoppv#1552)
  Workaround for slow scipy truncnorm by using the version from 1.3.3 (metoppv#1576)
  Support for a static additional predictor within the EMOS plugins (metoppv#1564)
  Fix negative grid spacing (metoppv#1583)
  Increase leniency of EMOS application (metoppv#1577)

Author: gavinevans

CONTRIBUTING.md

@@ -55,6 +55,7 @@ below:
 * Daniel Mentiplay (Bureau of Meteorology, Australia)
 * Stephen Moseley (Met Office, UK)
 * Meabh NicGuidhir (Met Office, UK)
+* Carwyn Pelley (Met Office, UK)
 * Tim Pillinger (Met Office, UK)
 * Fiona Rust (Met Office, UK)
 * Chris Sampson (Met Office, UK)

improver/calibration/ensemble_calibration.py

@@ -285,12 +285,29 @@ def merge_land_and_sea(calibrated_land_only: Cube, uncalibrated: Cube) -> None:
         calibrated_land_only.data = new_data
 
 
+def _ceiling_fp(cube: Cube) -> np.ndarray:
+    """Find the forecast period points rounded up to the next hour.
+
+    Args:
+        cube:
+            Cube with a forecast_period coordinate.
+
+    Returns:
+        The forecast period points in units of hours after
+        rounding the points up to the next hour.
+    """
+    coord = cube.coord("forecast_period").copy()
+    coord.convert_units("hours")
+    return np.ceil(coord.points)
+
+
 def forecast_coords_match(first_cube: Cube, second_cube: Cube) -> None:
     """
-    Determine if two cubes have equivalent forecast_periods and that the hours
-    of the forecast_reference_time coordinates match. Only the point of the
-    forecast reference time coordinate is checked to ensure that a calibration
-    / coefficient cube matches the forecast cube, as appropriate.
+    Determine if two cubes have equivalent forecast_periods and
+    forecast_reference_time coordinates with an accepted leniency.
+    The forecast period is rounded up to the next hour to
+    support calibrating subhourly forecasts with coefficients taken from on
+    the hour. For forecast reference time, only the hour is checked.
 
     Args:
         first_cube:
@@ -302,8 +319,8 @@ def forecast_coords_match(first_cube: Cube, second_cube: Cube) -> None:
         ValueError: The two cubes are not equivalent.
     """
     mismatches = []
-    if first_cube.coord("forecast_period") != second_cube.coord("forecast_period"):
-        mismatches.append("forecast_period")
+    if _ceiling_fp(first_cube) != _ceiling_fp(second_cube):
+        mismatches.append("rounded forecast_period hours")
 
     if get_frt_hours(first_cube.coord("forecast_reference_time")) != get_frt_hours(
         second_cube.coord("forecast_reference_time")

improver/calibration/utilities.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/cli/apply_emos_coefficients.py

@@ -75,8 +75,6 @@ def process(
     from improver.wxcode.weather_symbols import WeatherSymbols
 
     if not cubes:
-        raise RuntimeError(
-            "Not enough input arguments. " "See help for more information."
-        )
+        raise RuntimeError("Not enough input arguments. See help for more information.")
 
     return WeatherSymbols(wxtree, model_id_attr=model_id_attr)(CubeList(cubes))

improver/cli/wxcode.py

@@ -0,0 +1,59 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# -----------------------------------------------------------------------------
+# (C) British Crown Copyright 2017-2021 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 generate modal weather symbols over periods."""
+
+from improver import cli
+
+
+@cli.clizefy
+@cli.with_output
+def process(*cubes: cli.inputcube):
+    """Generates a modal weather symbol for the period covered by the input
+    weather symbol cubes. Where there are different weather codes available
+    for night and day, the modal code returned is always a day code, regardless
+    of the times covered by the input files.
+
+    Args:
+        cubes (iris.cube.CubeList):
+            A cubelist containing weather symbols cubes that cover the period
+            over which a modal symbol is desired.
+
+    Returns:
+        iris.cube.Cube:
+            A cube of modal weather symbols over a period.
+    """
+    from improver.wxcode.modal_code import ModalWeatherCode
+
+    if not cubes:
+        raise RuntimeError("Not enough input arguments. See help for more information.")
+
+    return ModalWeatherCode()(cubes)

improver/cli/wxcode_modal.py

@@ -0,0 +1,182 @@
+# -*- coding: utf-8 -*-
+# -----------------------------------------------------------------------------
+# (C) British Crown Copyright 2017-2021 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.
+"""
+This module defines the truncnorm as per scipy v1.3.3 to overcome performance
+issue introduced in later versions:
+- https://github.com/scipy/scipy/issues/12370
+- https://github.com/scipy/scipy/issues/12733
+
+"""
+import numpy as np
+import scipy.special as sc
+from scipy.stats._distn_infrastructure import rv_continuous
+
+# ============================================================================
+# |                        Copyright SciPy                                   |
+# | Code from this point unto the termination banner is copyright SciPy.     |
+# |                                                                          |
+# | Copyright © 2001, 2002 Enthought, Inc.                                   |
+# | All rights reserved.                                                     |
+# |                                                                          |
+# | Copyright © 2003-2019 SciPy Developers.                                  |
+# | 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 Enthought nor the names of the SciPy Developers    |
+# |   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 REGENTS 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.             |
+# |                                                                          |
+# | Further details can be found at scipy.org/scipylib/license.html          |
+# ============================================================================
+
+# Source: https://github.com/scipy/scipy/blob/v1.3.3/scipy/stats/_continuous_\
+# distns.py
+
+
+_norm_pdf_C = np.sqrt(2 * np.pi)
+_norm_pdf_logC = np.log(_norm_pdf_C)
+
+
+def _norm_pdf(x):
+    return np.exp(-(x ** 2) / 2.0) / _norm_pdf_C
+
+
+def _norm_logpdf(x):
+    return -(x ** 2) / 2.0 - _norm_pdf_logC
+
+
+def _norm_cdf(x):
+    return sc.ndtr(x)
+
+
+def _norm_ppf(q):
+    return sc.ndtri(q)
+
+
+def _norm_sf(x):
+    return _norm_cdf(-x)
+
+
+def _norm_isf(q):
+    return -_norm_ppf(q)
+
+
+class truncnorm_gen(rv_continuous):
+    r"""A truncated normal continuous random variable.
+
+    %(before_notes)s
+
+    Notes
+    -----
+    The standard form of this distribution is a standard normal truncated to
+    the range [a, b] --- notice that a and b are defined over the domain of the
+    standard normal.  To convert clip values for a specific mean and standard
+    deviation, use::
+
+        a, b = (myclip_a - my_mean) / my_std, (myclip_b - my_mean) / my_std
+
+    `truncnorm` takes :math:`a` and :math:`b` as shape parameters.
+
+    %(after_notes)s
+
+    %(example)s
+
+    """
+
+    def _argcheck(self, a, b):
+        return a < b
+
+    def _get_support(self, a, b):
+        return a, b
+
+    def _get_norms(self, a, b):
+        _nb = _norm_cdf(b)
+        _na = _norm_cdf(a)
+        _sb = _norm_sf(b)
+        _sa = _norm_sf(a)
+        _delta = np.where(a > 0, _sa - _sb, _nb - _na)
+        with np.errstate(divide="ignore"):
+            return _na, _nb, _sa, _sb, _delta, np.log(_delta)
+
+    def _pdf(self, x, a, b):
+        ans = self._get_norms(a, b)
+        _delta = ans[4]
+        return _norm_pdf(x) / _delta
+
+    def _logpdf(self, x, a, b):
+        ans = self._get_norms(a, b)
+        _logdelta = ans[5]
+        return _norm_logpdf(x) - _logdelta
+
+    def _cdf(self, x, a, b):
+        ans = self._get_norms(a, b)
+        _na, _delta = ans[0], ans[4]
+        return (_norm_cdf(x) - _na) / _delta
+
+    def _ppf(self, q, a, b):
+        # XXX Use _lazywhere...
+        ans = self._get_norms(a, b)
+        _na, _nb, _sa, _sb = ans[:4]
+        ppf = np.where(
+            a > 0,
+            _norm_isf(q * _sb + _sa * (1.0 - q)),
+            _norm_ppf(q * _nb + _na * (1.0 - q)),
+        )
+        return ppf
+
+
+truncnorm = truncnorm_gen(name="truncnorm")
+
+
+# ============================================================================
+# |                        END SciPy copyright                               |
+# ============================================================================

improver/ensemble_copula_coupling/_scipy_continuous_distns.py

@@ -42,6 +42,7 @@
 from numpy import ndarray
 from scipy import stats
 
+import improver.ensemble_copula_coupling._scipy_continuous_distns as scipy_cont_distns
 from improver import BasePlugin
 from improver.calibration.utilities import convert_cube_data_to_2d
 from improver.ensemble_copula_coupling.utilities import (
@@ -727,14 +728,18 @@ def __init__(
                 a lower bound of zero should be [0, np.inf].
 
         """
-        try:
-            self.distribution = getattr(stats, distribution)
-        except AttributeError as err:
-            msg = (
-                "The distribution requested {} is not a valid distribution "
-                "in scipy.stats. {}".format(distribution, err)
-            )
-            raise AttributeError(msg)
+        if distribution == "truncnorm":
+            # Use scipy v1.3.3 truncnorm
+            self.distribution = scipy_cont_distns.truncnorm
+        else:
+            try:
+                self.distribution = getattr(stats, distribution)
+            except AttributeError as err:
+                msg = (
+                    "The distribution requested {} is not a valid distribution "
+                    "in scipy.stats. {}".format(distribution, err)
+                )
+                raise AttributeError(msg)
 
         if shape_parameters is None:
             if self.distribution.name == "truncnorm":

improver/ensemble_copula_coupling/ensemble_copula_coupling.py

@@ -69,7 +69,9 @@ def calculate_input_grid_spacing(cube_in: Cube) -> Tuple[float, float]:
     lon_spacing = calculate_grid_spacing(cube_in, "degree", axis="x", rtol=4.0e-5)
     lat_spacing = calculate_grid_spacing(cube_in, "degree", axis="y", rtol=4.0e-5)
 
-    if lon_spacing < 0 or lat_spacing < 0:
+    y_coord = cube_in.coord(axis="y").points
+    x_coord = cube_in.coord(axis="x").points
+    if x_coord[-1] < x_coord[0] or y_coord[-1] < y_coord[0]:
         raise ValueError("Input grid coordinates are not ascending.")
     return lat_spacing, lon_spacing
 

improver/regrid/grid.py

@@ -80,7 +80,8 @@ def calculate_grid_spacing(
     cube: Cube, units: Union[Unit, str], axis: str = "x", rtol: float = 1.0e-5
 ) -> float:
     """
-    Returns the grid spacing of a given spatial axis
+    Returns the grid spacing of a given spatial axis. This will be positive for
+    axes that stride negatively.
 
     Args:
         cube:
@@ -100,7 +101,7 @@ def calculate_grid_spacing(
     """
     coord = cube.coord(axis=axis).copy()
     coord.convert_units(units)
-    diffs = np.diff(coord.points)
+    diffs = np.abs(np.diff(coord.points))
     diffs_mean = np.mean(diffs)
 
     if not np.allclose(diffs, diffs_mean, rtol=rtol, atol=0.0):