Repository: pySTEPS/pysteps (language: python)

Issue:
Implement linear blending
Implement a linear blending algorithm which can be applied on the different nowcast algorithms.



Discussion / hints:

Fix patch:
diff --git a/doc/source/pysteps_reference/blending.rst b/doc/source/pysteps_reference/blending.rst
new file mode 100644
index 000000000..b827e1aa7
--- /dev/null
+++ b/doc/source/pysteps_reference/blending.rst
@@ -0,0 +1,8 @@
+================
+pysteps.blending
+================
+
+Implementation of blending methods for blending (ensemble) nowcasts with Numerical Weather Prediction (NWP) models.
+
+
+.. automodule:: pysteps.blending.linear_blending
diff --git a/examples/plot_linear_blending.py b/examples/plot_linear_blending.py
new file mode 100644
index 000000000..dada60256
--- /dev/null
+++ b/examples/plot_linear_blending.py
@@ -0,0 +1,190 @@
+# -*- coding: utf-8 -*-
+
+"""
+Plot linear blending
+====================
+
+This tutorial shows how to construct a simple linear blending between an ensemble 
+nowcast and a dummy Numerical Weather Prediction (NWP) rainfall forecast.
+"""
+
+# TODO: Make xarray ready for pysteps-v2
+
+from matplotlib import pyplot as plt
+import numpy as np
+from pprint import pprint
+from datetime import datetime
+from pysteps.motion.lucaskanade import dense_lucaskanade
+from pysteps import io, rcparams
+from pysteps.utils import conversion, transformation
+from pysteps.visualization import plot_precip_field
+from pysteps.utils import dimension
+from pysteps.blending.linear_blending import forecast
+
+
+def gaussian(x, max, mean, sigma):
+    return max * np.exp(-(x - mean) * (x - mean) / sigma / sigma / 2)
+
+
+def dummy_nwp(precip, n_leadtimes, max=20, mean=0, sigma=0.25, speed=100):
+    """Generates dummy NWP data with the same dimension as the input
+    precipitation field precip. The NWP data is a vertical line with a Gaussian
+    profile moving to the left"""
+
+    # precip is original radar image
+    rows = precip.shape[0]
+    cols = precip.shape[1]
+
+    # Initialise the dummy NWP data
+    precip_nwp = np.zeros((n_leadtimes, rows, cols))
+    x = np.linspace(-5, 5, cols)
+
+    for n in range(n_leadtimes):
+        for i in range(rows):
+            precip_nwp[n, i, :] = gaussian(x, max, mean, sigma)
+        mean -= speed / rows
+
+    return precip_nwp
+
+
+###############################################################################
+# Set nowcast parameters
+n_ens_members = 10
+n_leadtimes = 18
+start_blending = 30  # in minutes
+end_blending = 60  # in minutes
+seed = 24
+
+# Read precipitation field
+# ------------------------
+#
+# First thing, the sequence of Swiss radar composites is imported, converted and
+# transformed into units of dBR.
+
+date = datetime.strptime("201701311200", "%Y%m%d%H%M")
+data_source = "mch"
+
+# Load data source config
+root_path = rcparams.data_sources[data_source]["root_path"]
+path_fmt = rcparams.data_sources[data_source]["path_fmt"]
+fn_pattern = rcparams.data_sources[data_source]["fn_pattern"]
+fn_ext = rcparams.data_sources[data_source]["fn_ext"]
+importer_name = rcparams.data_sources[data_source]["importer"]
+importer_kwargs = rcparams.data_sources[data_source]["importer_kwargs"]
+timestep = rcparams.data_sources[data_source]["timestep"]
+
+# Find the radar files in the archive
+fns = io.find_by_date(
+    date, root_path, path_fmt, fn_pattern, fn_ext, timestep, num_prev_files=2
+)
+
+# Read the data from the archive
+importer = io.get_method(importer_name, "importer")
+precip = io.read_timeseries(fns, importer, legacy=False, **importer_kwargs)
+
+# Get the metadata
+metadata = precip.x.attrs.copy()
+metadata.update(**precip.y.attrs)
+metadata.update(**precip.attrs)
+
+# Convert to rain rate
+precip, metadata = conversion.to_rainrate(precip, metadata)
+
+# Upscale data to 2 km to limit memory usage
+precip, metadata = dimension.aggregate_fields_space(precip, metadata, 2000)
+
+# Import the dummy NWP data (vertical line moving to the left)
+precip_nwp = dummy_nwp(
+    precip[-1, :, :], n_leadtimes + 1, max=7, mean=4, speed=0.2 * 350
+)
+metadata_nwp = metadata.copy()
+
+# Plot the radar rainfall field and the first time step of the dummy NWP forecast.
+date_str = datetime.strftime(date, "%Y-%m-%d %H:%M")
+plt.figure(figsize=(10, 5))
+plt.subplot(121)
+plot_precip_field(
+    precip[-1, :, :], geodata=metadata, title=f"Radar observation at {date_str}"
+)
+plt.subplot(122)
+plot_precip_field(
+    precip_nwp[0, :, :], geodata=metadata_nwp, title=f"Dummy NWP forecast at {date_str}"
+)
+plt.tight_layout()
+plt.show()
+
+# Log-transform the data to unit of dBR, set the threshold to 0.1 mm/h,
+# set the fill value to -15 dBR
+precip, metadata = transformation.dB_transform(
+    precip, metadata, threshold=0.1, zerovalue=-15.0
+)
+
+# Nicely print the metadata
+pprint(metadata)
+
+###############################################################################
+# Linear blending of nowcast and NWP data
+
+# Estimate the motion field
+velocity = dense_lucaskanade(precip)
+
+# Define nowcast keyword arguments
+nowcast_kwargs = {
+    "n_ens_members": n_ens_members,
+    "n_cascade_levels": 6,
+    "R_thr": -10.0,
+    "kmperpixel": 2.0,
+    "timestep": timestep,
+    "noise_method": "nonparametric",
+    "vel_pert_method": "bps",
+    "mask_method": "incremental",
+}
+
+# Calculate the blended precipitation field
+precip_blended = forecast(
+    precip=precip[-3:, :, :],
+    precip_metadata=metadata,
+    velocity=velocity,
+    timesteps=n_leadtimes,
+    timestep=timestep,
+    nowcast_method="steps",
+    precip_nwp=precip_nwp[1:, :, :],
+    precip_nwp_metadata=metadata_nwp,
+    start_blending=start_blending,
+    end_blending=end_blending,
+    nowcast_kwargs=nowcast_kwargs,
+)
+
+"""
+extrap_kwargs = {"allow_nonfinite_values": True}
+nowcast_kwargs = {"extrap_kwargs": extrap_kwargs}
+
+# Calculate the blended precipitation field
+precip_blended = forecast(
+    precip[-1, :, :],
+    V,
+    n_leadtimes,
+    timestep,
+    "extrapolation",
+    precip_nwp=precip_nwp[1: :, :],
+    start_blending=start_blending,
+    end_blending=end_blending,
+    nowcast_kwargs=nowcast_kwargs,
+)
+"""
+
+# Calculate the ensemble average
+if len(precip_blended.shape) == 4:
+    precip_blended_mean = np.mean(precip_blended[:, :, :, :], axis=0)
+else:
+    precip_blended_mean = np.copy(precip_blended)
+
+# Plot the blended field
+for i in range(0, n_leadtimes, 3):
+    plot_precip_field(
+        precip_blended_mean[i, :, :],
+        geodata=metadata,
+        title="Blended field (+ %i min)" % ((i + 1) * timestep),
+    )
+    plt.show()
+    plt.close()
diff --git a/pysteps/blending/__init__.py b/pysteps/blending/__init__.py
new file mode 100644
index 000000000..420a11bd0
--- /dev/null
+++ b/pysteps/blending/__init__.py
@@ -0,0 +1,4 @@
+# -*- coding: utf-8 -*-
+"""Methods for blending NWP model(s) with nowcasts."""
+
+from pysteps.blending.interface import get_method
diff --git a/pysteps/blending/interface.py b/pysteps/blending/interface.py
new file mode 100644
index 000000000..6764e95ba
--- /dev/null
+++ b/pysteps/blending/interface.py
@@ -0,0 +1,48 @@
+# -*- coding: utf-8 -*-
+"""
+pysteps.blending.interface
+==========================
+Interface for the blending module. It returns a callable function for computing
+blended nowcasts with NWP models.
+.. autosummary::
+    :toctree: ../generated/
+    get_method
+"""
+
+from pysteps.blending import linear_blending
+
+_blending_methods = dict()
+_blending_methods["linear_blending"] = linear_blending.forecast
+
+
+def get_method(name):
+    """Return a callable function for computing nowcasts.
+    Description:
+    Return a callable function for computing deterministic or ensemble
+    precipitation nowcasts.
+    Implemented methods:
+    +-----------------+-------------------------------------------------------+
+    |     Name        |              Description                              |
+    +=================+=======================================================+
+    +-----------------+-------------------------------------------------------+
+    |  linear         | the linear blending of a nowcast method with other    |
+    |  blending       | data (e.g. NWP data).                                 |
+    +-----------------+-------------------------------------------------------+
+    """
+    if isinstance(name, str):
+        name = name.lower()
+    else:
+        raise TypeError(
+            "Only strings supported for the method's names.\n"
+            + "Available names:"
+            + str(list(_blending_methods.keys()))
+        ) from None
+
+    try:
+        return _blending_methods[name]
+    except KeyError:
+        raise ValueError(
+            "Unknown blending method {}\n".format(name)
+            + "The available methods are:"
+            + str(list(_blending_methods.keys()))
+        ) from None
diff --git a/pysteps/blending/linear_blending.py b/pysteps/blending/linear_blending.py
new file mode 100644
index 000000000..5b8ea5a21
--- /dev/null
+++ b/pysteps/blending/linear_blending.py
@@ -0,0 +1,200 @@
+# -*- coding: utf-8 -*-
+"""
+pysteps.nowcasts.linear_blending
+================================
+
+Implementation of the linear blending between nowcast and NWP data.
+
+.. autosummary::
+    :toctree: ../generated/
+
+    forecast
+"""
+
+import numpy as np
+from pysteps import nowcasts
+from pysteps.utils import conversion
+
+# TODO: Make xarray ready for pysteps-v2
+
+
+def forecast(
+    precip,
+    precip_metadata,
+    velocity,
+    timesteps,
+    timestep,
+    nowcast_method,
+    precip_nwp=None,
+    precip_nwp_metadata=None,
+    start_blending=120,
+    end_blending=240,
+    fill_nwp=True,
+    nowcast_kwargs=dict(),
+):
+
+    """Generate a forecast by linearly blending nowcasts with NWP data
+
+    Parameters
+    ----------
+    precip: array_like
+      Array containing the input precipitation field(s) ordered by timestamp
+      from oldest to newest. The time steps between the inputs are assumed
+      to be regular.
+    precip_metadata: dict
+        Metadata dictionary containing (at least) the transform, unit and threshold
+        attributes as described in the documentation of :py:mod:`pysteps.io.importers`.
+    velocity; array_like
+      Array of shape (2, m, n) containing the x- and y-components of the advection
+      field. The velocities are assumed to represent one time step between the
+      inputs. All values are required to be finite.
+    timesteps: int
+      Number of time steps to forecast.
+    timestep: int or float
+      The time difference (in minutes) between consecutive forecast fields.
+    nowcast_method: {'anvil', 'eulerian', 'lagrangian', 'extrapolation', 'lagrangian_probability', 'sprog', 'steps', 'sseps'}
+      Name of the nowcast method used to forecast the precipitation.
+    precip_nwp: array_like or NoneType, optional
+      Array of shape (timesteps, m, n) in the case of no ensemble or
+      of shape (n_ens_members, timesteps, m, n) in the case of an ensemble
+      containing the NWP precipitation fields ordered by timestamp from oldest
+      to newest. The time steps between the inputs are assumed to be regular
+      (and identical to the time step between the nowcasts). If no NWP
+      data is given the value of precip_nwp is None and no blending will be performed.
+    precip_nwp_metadata: dict or NoneType, optional
+        NWP metadata dictionary containing (at least) the transform, unit and threshold
+        attributes as described in the documentation of :py:mod:`pysteps.io.importers`.
+    start_blending: int, optional
+      Time stamp (in minutes) after which the blending should start. Before this
+      only the nowcast data is used.
+    end_blending: int, optional
+      Time stamp (in minutes) after which the blending should end. Between
+      start_blending and end_blending the nowcasts and NWP data are linearly
+      merged with each other. After end_blending only the NWP data is used.
+    fill_nwp: bool, optional
+      Standard value is True. If True, the NWP data will be used to fill in the
+      no data mask of the nowcast.
+    nowcast_kwargs: dict, optional
+      Dictionary containing keyword arguments for the nowcast method.
+
+    Returns
+    -------
+    R_blended: ndarray
+      Array of shape (timesteps, m, n) in the case of no ensemble or
+      of shape (n_ens_members, timesteps, m, n) in the case of an ensemble
+      containing the precipation forecast generated by linearly blending
+      the nowcasts and the NWP data. n_ens_members equals the maximum no. of
+      ensemble members in either the nowcast or nwp model(s).
+    """
+
+    # Calculate the nowcasts
+    nowcast_method_func = nowcasts.get_method(nowcast_method)
+    precip_nowcast = nowcast_method_func(
+        precip,
+        velocity,
+        timesteps,
+        **nowcast_kwargs,
+    )
+
+    # Make sure that precip_nowcast and precip_nwp are in mm/h
+    precip_nowcast, _ = conversion.to_rainrate(precip_nowcast, metadata=precip_metadata)
+    precip_nwp, _ = conversion.to_rainrate(precip_nwp, metadata=precip_nwp_metadata)
+
+    # Check if NWP data is given as input
+    if precip_nwp is not None:
+
+        if len(precip_nowcast.shape) == 4:
+            n_ens_members_nowcast = precip_nowcast.shape[0]
+            if n_ens_members_nowcast == 1:
+                precip_nowcast = np.squeeze(precip_nowcast)
+        else:
+            n_ens_members_nowcast = 1
+
+        if len(precip_nwp.shape) == 4:
+            n_ens_members_nwp = precip_nwp.shape[0]
+            if n_ens_members_nwp == 1:
+                precip_nwp = np.squeeze(precip_nwp)
+        else:
+            n_ens_members_nwp = 1
+
+        # Now, repeat the nowcast ensemble members or the nwp models/members until
+        # it has the same amount of members as n_ens_members_max. For instance, if
+        # you have 10 ensemble nowcasts members and 3 NWP members, the output will
+        # be an ensemble of 10 members. Hence, the three NWP members are blended
+        # with the first three members of the nowcast (member one with member one,
+        # two with two, etc.), subsequently, the same NWP members are blended with
+        # the next three members (NWP member one with member 4, NWP member 2 with
+        # member 5, etc.), until 10 is reached.
+        n_ens_members_max = max(n_ens_members_nowcast, n_ens_members_nwp)
+        n_ens_members_min = min(n_ens_members_nowcast, n_ens_members_nwp)
+
+        if n_ens_members_min != n_ens_members_max:
+            if n_ens_members_nwp == 1:
+                precip_nwp = np.repeat(
+                    precip_nwp[np.newaxis, :, :], n_ens_members_max, axis=0
+                )
+            elif n_ens_members_nowcast == 1:
+                precip_nowcast = np.repeat(
+                    precip_nowcast[np.newaxis, :, :], n_ens_members_max, axis=0
+                )
+            else:
+                repeats = [
+                    (n_ens_members_max + i) // n_ens_members_min
+                    for i in range(n_ens_members_min)
+                ]
+
+                if n_ens_members_nwp == n_ens_members_min:
+                    precip_nwp = np.repeat(precip_nwp, repeats, axis=0)
+                elif n_ens_members_nowcast == n_ens_members_min:
+                    precip_nowcast = np.repeat(precip_nowcast, repeats, axis=0)
+
+        # Check if dimensions are correct
+        assert (
+            precip_nwp.shape == precip_nowcast.shape
+        ), "The dimensions of precip_nowcast and precip_nwp need to be identical: dimension of precip_nwp = {} and dimension of precip_nowcast = {}".format(
+            precip_nwp.shape, precip_nowcast.shape
+        )
+
+        # Initialise output
+        R_blended = np.zeros_like(precip_nowcast)
+
+        # Calculate the weights
+        for i in range(timesteps):
+            # Calculate what time we are at
+            t = (i + 1) * timestep
+
+            # Calculate the weight with a linear relation (weight_nwp at start_blending = 0.0)
+            # and (weight_nwp at end_blending = 1.0)
+            weight_nwp = (t - start_blending) / (end_blending - start_blending)
+
+            # Set weights at times before start_blending and after end_blending
+            if weight_nwp < 0.0:
+                weight_nwp = 0.0
+            elif weight_nwp > 1.0:
+                weight_nwp = 1.0
+
+            # Calculate weight_nowcast
+            weight_nowcast = 1.0 - weight_nwp
+
+            # Calculate output by combining precip_nwp and precip_nowcast,
+            # while distinguishing between ensemble and non-ensemble methods
+            if n_ens_members_max == 1:
+                R_blended[i, :, :] = (
+                    weight_nwp * precip_nwp[i, :, :]
+                    + weight_nowcast * precip_nowcast[i, :, :]
+                )
+            else:
+                R_blended[:, i, :, :] = (
+                    weight_nwp * precip_nwp[:, i, :, :]
+                    + weight_nowcast * precip_nowcast[:, i, :, :]
+                )
+
+            # Find where the NaN values are and replace them with NWP data
+            if fill_nwp:
+                nan_indices = np.isnan(R_blended)
+                R_blended[nan_indices] = precip_nwp[nan_indices]
+    else:
+        # If no NWP data is given, the blended field is simply equal to the nowcast field
+        R_blended = precip_nowcast
+
+    return R_blended
diff --git a/pysteps/nowcasts/interface.py b/pysteps/nowcasts/interface.py
index a86e6c080..0b7b94bce 100644
--- a/pysteps/nowcasts/interface.py
+++ b/pysteps/nowcasts/interface.py
@@ -32,7 +32,14 @@
 """
 
 from pysteps.extrapolation.interface import eulerian_persistence
-from pysteps.nowcasts import anvil, extrapolation, linda, sprog, steps, sseps
+from pysteps.nowcasts import (
+    anvil,
+    extrapolation,
+    linda,
+    sprog,
+    steps,
+    sseps,
+)
 from pysteps.nowcasts import lagrangian_probability
 
 _nowcast_methods = dict()
diff --git a/pysteps/nowcasts/lagrangian_probability.py b/pysteps/nowcasts/lagrangian_probability.py
index c9889da6b..2bc8a5a08 100644
--- a/pysteps/nowcasts/lagrangian_probability.py
+++ b/pysteps/nowcasts/lagrangian_probability.py
@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
 """
 pysteps.nowcasts.lagrangian_probability
-======================================
+=======================================
 
 Implementation of the local Lagrangian probability nowcasting technique
 described in :cite:`GZ2004`.


Test patch:
diff --git a/pysteps/tests/test_blending_linear_blending.py b/pysteps/tests/test_blending_linear_blending.py
new file mode 100644
index 000000000..712963eba
--- /dev/null
+++ b/pysteps/tests/test_blending_linear_blending.py
@@ -0,0 +1,109 @@
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import pytest
+from pysteps.blending.linear_blending import forecast
+from numpy.testing import assert_array_almost_equal
+from pysteps.utils import transformation
+
+# Test function arguments
+linear_arg_values = [
+    (5, 30, 60, 20, 45, "eulerian", None),
+    (4, 23, 33, 9, 28, "eulerian", None),
+    (3, 18, 36, 13, 27, "eulerian", None),
+    (7, 30, 68, 11, 49, "eulerian", None),
+    (10, 100, 160, 25, 130, "eulerian", None),
+    (6, 60, 180, 22, 120, "eulerian", None),
+    (5, 100, 200, 40, 150, "eulerian", None),
+    (5, 30, 60, 20, 45, "extrapolation", np.zeros((2, 200, 200))),
+    (4, 23, 33, 9, 28, "extrapolation", np.zeros((2, 200, 200))),
+    (3, 18, 36, 13, 27, "extrapolation", np.zeros((2, 200, 200))),
+    (7, 30, 68, 11, 49, "extrapolation", np.zeros((2, 200, 200))),
+    (10, 100, 160, 25, 130, "extrapolation", np.zeros((2, 200, 200))),
+    (6, 60, 180, 22, 120, "extrapolation", np.zeros((2, 200, 200))),
+    (5, 100, 200, 40, 150, "extrapolation", np.zeros((2, 200, 200))),
+]
+
+
+@pytest.mark.parametrize(
+    "timestep, start_blending, end_blending, n_timesteps, controltime, nowcast_method, V",
+    linear_arg_values,
+)
+def test_linear_blending(
+    timestep, start_blending, end_blending, n_timesteps, controltime, nowcast_method, V
+):
+    """Tests if the linear blending function is correct. For the nowcast data a precipitation field
+    which is constant over time is taken. One half of the field has no rain and the other half
+    has a set value. For the NWP data a similar field is taken, the only difference
+    being that now the other half of the field is zero. The blended field should have a
+    constant value over the entire field at the timestep right in the middle between the start
+    of the blending and the end of the blending. This assertion is checked to see if the
+    linear blending function works well."""
+
+    # The argument controltime gives the timestep at which the field is assumed to be
+    # entirely constant
+
+    # Assert that the control time step is in the range of the forecasted time steps
+    assert controltime <= (
+        n_timesteps * timestep
+    ), "Control time needs to be within reach of forecasts, controltime = {} and n_timesteps = {}".format(
+        controltime, n_timesteps
+    )
+
+    # Assert that the start time of the blending comes before the end time of the blending
+    assert (
+        start_blending < end_blending
+    ), "Start time of blending needs to be smaller than end time of blending"
+
+    # Assert that the control time is a multiple of the time step
+    assert (
+        not controltime % timestep
+    ), "Control time needs to be a multiple of the time step"
+
+    # Initialise dummy NWP data
+    r_nwp = np.zeros((n_timesteps, 200, 200))
+
+    for i in range(100):
+        r_nwp[:, i, :] = 11.0
+
+    # Define nowcast input data
+    r_input = np.zeros((200, 200))
+
+    for i in range(100, 200):
+        r_input[i, :] = 11.0
+
+    # Transform from mm/h to dB
+    r_input, _ = transformation.dB_transform(
+        r_input, None, threshold=0.1, zerovalue=-15.0
+    )
+
+    # Calculate the blended field
+    r_blended = forecast(
+        r_input,
+        dict({"unit": "mm/h", "transform": "dB"}),
+        V,
+        n_timesteps,
+        timestep,
+        nowcast_method,
+        r_nwp,
+        dict({"unit": "mm/h", "transform": None}),
+        start_blending=start_blending,
+        end_blending=end_blending,
+    )
+
+    # Assert that the blended field has the expected dimension
+    assert r_blended.shape == (
+        n_timesteps,
+        200,
+        200,
+    ), "The shape of the blended array does not have the expected value. The shape is {}".format(
+        r_blended.shape
+    )
+
+    # Assert that the blended field at the control time step is equal to
+    # a constant field with the expected value.
+    assert_array_almost_equal(
+        r_blended[controltime // timestep - 1],
+        np.ones((200, 200)) * 5.5,
+        err_msg="The blended array does not have the expected value",
+    )