Linear blending

doc/source/pysteps_reference/nowcasts.rst

@@ -14,3 +14,4 @@ Implementation of deterministic and ensemble nowcasting methods.
 .. automodule:: pysteps.nowcasts.sseps
 .. automodule:: pysteps.nowcasts.steps
 .. automodule:: pysteps.nowcasts.utils
+.. automodule:: pysteps.nowcasts.linear_blending

examples/plot_linear_blending.py

@@ -0,0 +1,162 @@
+# -*- coding: utf-8 -*-
+
+"""
+Plot linear blending
+====================
+"""
+
+from matplotlib import cm, pyplot as plt
+import numpy as np
+from pprint import pprint
+from pysteps.motion.lucaskanade import dense_lucaskanade
+from pysteps import io, rcparams, nowcasts
+from pysteps.utils import conversion, transformation
+from pysteps.visualization import plot_precip_field
+from datetime import datetime
+from pysteps.utils import dimension
+from pysteps.nowcasts.linear_blending import forecast
+
+
+def gaussian(x, max, mean, sigma):
+    return max * np.exp(-(x - mean) * (x - mean) / sigma / sigma / 2)
+
+
+def dummy_nwp(R, n_leadtimes, max=20, mean=0, sigma=0.25, speed=100):
+    """Generates dummy NWP data with the same dimension as the input
+    precipitation field R. The NWP data is a vertical line with a Gaussian
+    profile moving to the left"""
+
+    # R is original radar image
+    rows = R.shape[0]
+    cols = R.shape[1]
+
+    # Initialise the dummy NWP data
+    R_nwp = np.zeros((n_leadtimes, rows, cols))
+    x = np.linspace(-5, 5, cols)
+
+    for n in range(n_leadtimes):
+        for i in range(rows):
+            R_nwp[n, i, :] = gaussian(x, max, mean, sigma)
+        mean -= speed / rows
+
+    return R_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")
+R, _, metadata = io.read_timeseries(fns, importer, legacy=True, **importer_kwargs)
+
+# Convert to rain rate
+R, metadata = conversion.to_rainrate(R, metadata)
+
+# Upscale data to 2 km to limit memory usage
+R, metadata = dimension.aggregate_fields_space(R, metadata, 2000)
+
+# Import the dummy NWP data (vertical line moving to the left)
+R_nwp = dummy_nwp(R[-1, :, :], n_leadtimes + 1, max=7, mean=4, speed=0.2 * 350)
+
+# Log-transform the data to unit of dBR, set the threshold to 0.1 mm/h,
+# set the fill value to -15 dBR
+R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0)
+
+# Set missing values with the fill value
+# R[~np.isfinite(R)] = -15.0
+
+# Nicely print the metadata
+pprint(metadata)
+
+###############################################################################
+# Linear blending of nowcast and NWP data
+
+# Estimate the motion field
+V = dense_lucaskanade(R)
+
+# 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
+R_blended = forecast(
+    R[-3:, :, :],
+    V,
+    n_leadtimes,
+    timestep,
+    "steps",
+    R_nwp=R_nwp[1:, :, :],
+    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
+R_blended = forecast(
+    R[-1, :, :],
+    V,
+    n_leadtimes,
+    timestep,
+    "extrapolation",
+    R_nwp=R_nwp[1: :, :],
+    start_blending=start_blending,
+    end_blending=end_blending,
+    nowcast_kwargs=nowcast_kwargs,
+)
+"""
+
+# Calculate the ensemble average
+if len(R_blended.shape) == 4:
+    R_blended_mean = np.mean(R_blended[:, :, :, :], axis=0)
+else:
+    R_blended_mean = np.copy(R_blended)
+
+# Plot the blended field
+for i in range(n_leadtimes):
+    plot_precip_field(
+        R_blended_mean[i, :, :],
+        geodata=metadata,
+        title="Blended field (+ %i min)" % ((i + 1) * timestep),
+    )
+    plt.show()
+    plt.close()

pysteps/nowcasts/interface.py

@@ -32,7 +32,15 @@
 """
 
 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,
+    linear_blending,
+)
 from pysteps.nowcasts import lagrangian_probability
 
 _nowcast_methods = dict()
@@ -46,6 +54,7 @@
 _nowcast_methods["sprog"] = sprog.forecast
 _nowcast_methods["sseps"] = sseps.forecast
 _nowcast_methods["steps"] = steps.forecast
+_nowcast_methods["linear blending"] = linear_blending.forecast
 
 
 def get_method(name):
@@ -83,6 +92,9 @@ def get_method(name):
     |  sseps          | short-space ensemble prediction system (SSEPS).       |
     |                 | Essentially, this is a localization of STEPS          |
     +-----------------+-------------------------------------------------------+
+    |  linear         | the linear blending of a nowcast method with other    |
+    |  blending       | data (e.g. NWP data).                                 |
+    +-----------------+-------------------------------------------------------+
     """
     if isinstance(name, str):
         name = name.lower()

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`.

pysteps/nowcasts/linear_blending.py

@@ -0,0 +1,172 @@
+# -*- 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 transformation
+
+
+def forecast(
+    precip,
+    velocity,
+    timesteps,
+    timestep,
+    nowcast_method,
+    R_nwp=None,
+    start_blending=120,
+    end_blending=240,
+    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.
+    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.
+    R_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 R_nwp is None and no blending will be performed.
+    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.
+    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.
+    """
+
+    # Calculate the nowcasts
+    nowcast_method_func = nowcasts.get_method(nowcast_method)
+    R_nowcast = nowcast_method_func(
+        precip,
+        velocity,
+        timesteps,
+        **nowcast_kwargs,
+    )
+
+    # Transform the precipitation back to mm/h
+    R_nowcast = transformation.dB_transform(R_nowcast, threshold=-10.0, inverse=True)[0]
+
+    # Check if NWP data is given as input
+    if R_nwp is not None:
+
+        if len(R_nowcast.shape) == 4:
+            n_ens_members_nowcast = R_nowcast.shape[0]
+            if n_ens_members_nowcast == 1:
+                R_nowcast = np.squeeze(R_nowcast)
+        else:
+            n_ens_members_nowcast = 1
+
+        if len(R_nwp.shape) == 4:
+            n_ens_members_nwp = R_nwp.shape[0]
+            if n_ens_members_nwp == 1:
+                R_nwp = np.squeeze(R_nwp)
+        else:
+            n_ens_members_nwp = 1
+
+        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:
+                R_nwp = np.repeat(R_nwp[np.newaxis, :, :], n_ens_members_max, axis=0)
+            elif n_ens_members_nowcast == 1:
+                R_nowcast = np.repeat(
+                    R_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:
+                    R_nwp = np.repeat(R_nwp, repeats, axis=0)
+                elif n_ens_members_nowcast == n_ens_members_min:
+                    R_nowcast = np.repeat(R_nowcast, repeats, axis=0)
+
+        # Check if dimensions are correct
+        assert (
+            R_nwp.shape == R_nowcast.shape
+        ), "The dimensions of R_nowcast and R_nwp need to be identical: dimension of R_nwp = {} and dimension of R_nowcast = {}".format(
+            R_nwp.shape, R_nowcast.shape
+        )
+
+        # Initialise output
+        R_blended = np.zeros_like(R_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 R_nwp and R_nowcast,
+            # while distinguishing between ensemble and non-ensemble methods
+            if n_ens_members_max == 1:
+                R_blended[i, :, :] = (
+                    weight_nwp * R_nwp[i, :, :] + weight_nowcast * R_nowcast[i, :, :]
+                )
+            else:
+                R_blended[:, i, :, :] = (
+                    weight_nwp * R_nwp[:, i, :, :]
+                    + weight_nowcast * R_nowcast[:, i, :, :]
+                )
+
+            # Find where the NaN values are and replace them with NWP data
+            nan_indices = np.isnan(R_blended)
+            R_blended[nan_indices] = R_nwp[nan_indices]
+    else:
+        # If no NWP data is given, the blended field is simply equal to the nowcast field
+        R_blended = R_nowcast
+
+    return R_blended