Implement option for saliency-based blending (#285)

* Implement option for salience weight

* Remove code duplication and add context

* extend tests, add example and change docstring

* Improv example a

Co-authored-by: Nathalie Rombeek <nathalie.rombeek@meteoswiss.ch>

doc/source/references.bib

@@ -136,6 +136,17 @@ @ARTICLE{Her2000
   DOI = "10.1175/1520-0434(2000)015<0559:DOTCRP>2.0.CO;2"
 }
 
+@article{Hwang2015,
+  AUTHOR = "Hwang, Yunsung and Clark, Adam J and Lakshmanan, Valliappa and Koch, Steven E",
+  TITLE = "Improved nowcasts by blending extrapolation and model forecasts",
+  JOURNAL = "Weather and Forecasting",
+  VOLUME = 30,
+  NUMBER = 5,
+  PAGES = "1201--1217",
+  YEAR = 2015,
+  DOI = "10.1175/WAF-D-15-0057.1"
+}
+
 @ARTICLE{NBSG2017,
   AUTHOR = "D. Nerini and N. Besic and I. Sideris and U. Germann and L. Foresti",
   TITLE = "A non-stationary stochastic ensemble generator for radar rainfall fields based on the short-space {F}ourier transform",

examples/plot_linear_blending.py

@@ -15,7 +15,8 @@
 from matplotlib import pyplot as plt
 
 import pysteps
-from pysteps import io, rcparams, blending
+from pysteps import io, rcparams, nowcasts, blending
+from pysteps.utils import conversion
 from pysteps.visualization import plot_precip_field
 
 
@@ -149,41 +150,110 @@
 
 
 ################################################################################
-# Visualize the output
-# ~~~~~~~~~~~~~~~~~~~~
+# The salient blending of nowcast and NWP rainfall forecast
+# ---------------------------------------------------------
 #
+# This method follows the saliency-based blending procedure described in :cite:`Hwang2015`. The
+# blending is based on intensities and forecast times. The blended product preserves pixel
+# intensities with time if they are strong enough based on their ranked salience. Saliency is
+# the property of an object to be outstanding with respect to its surroundings. The ranked salience
+# is calculated by first determining the difference in the normalized intensity of the nowcasts
+# and NWP. Next, the pixel intensities are ranked, in which equally comparable values receive
+# the same ranking number.
+
+# Calculate the salient blended precipitation field
+precip_salient_blended = blending.linear_blending.forecast(
+    precip=radar_precip[-1, :, :],
+    precip_metadata=radar_metadata,
+    velocity=velocity_radar,
+    timesteps=18,
+    timestep=10,
+    nowcast_method="extrapolation",  # simple advection nowcast
+    precip_nwp=precip_nwp,
+    precip_nwp_metadata=nwp_metadata,
+    start_blending=60,  # in minutes (this is an arbritrary choice)
+    end_blending=120,  # in minutes (this is an arbritrary choice)
+    saliency=True,
+)
+
+
+################################################################################
+# Visualize the output
+# --------------------
+
+################################################################################
+# Calculate the radar rainfall nowcasts for visualization
+
+nowcast_method_func = nowcasts.get_method("extrapolation")
+precip_nowcast = nowcast_method_func(
+    precip=radar_precip[-1, :, :],
+    velocity=velocity_radar,
+    timesteps=18,
+)
+
+# Make sure that precip_nowcast are in mm/h
+precip_nowcast, _ = conversion.to_rainrate(precip_nowcast, metadata=radar_metadata)
+
+################################################################################
 # The linear blending starts at 60 min, so during the first 60 minutes the
 # blended forecast only consists of the extrapolation forecast (consisting of an
 # extrapolation nowcast). Between 60 and 120 min, the NWP forecast gradually gets more
-# weight, whereas the extrapolation forecasts gradually gets less weight.
-# After 120 min, the blended forecast entirely consists of the NWP rainfall
-# forecast.
+# weight, whereas the extrapolation forecasts gradually gets less weight. In addition,
+# the saliency-based blending takes also the difference in pixel intensities into account,
+# which are preserved over time if they are strong enough based on their ranked salience.
+# Furthermore, pixels with relative low intensities get a lower weight and stay smaller in
+# the saliency-based blending compared to linear blending. After 120 min, the blended
+# forecast entirely consists of the NWP rainfall forecast.
 
-fig = plt.figure(figsize=(4, 8))
+fig = plt.figure(figsize=(8, 12))
 
-leadtimes_min = [30, 60, 90, 120]
+leadtimes_min = [30, 60, 80, 100, 120]
 n_leadtimes = len(leadtimes_min)
 for n, leadtime in enumerate(leadtimes_min):
 
+    # Extrapolation
+    plt.subplot(n_leadtimes, 4, n * 4 + 1)
+    plot_precip_field(
+        precip_nowcast[int(leadtime / timestep) - 1, :, :],
+        geodata=radar_metadata,
+        title=f"Nowcast + {leadtime} min",
+        axis="off",
+        colorbar=False,
+    )
+
     # Nowcast with blending into NWP
-    plt.subplot(n_leadtimes, 2, n * 2 + 1)
+    plt.subplot(n_leadtimes, 4, n * 4 + 2)
     plot_precip_field(
         precip_blended[int(leadtime / timestep) - 1, :, :],
         geodata=radar_metadata,
-        title=f"Nowcast +{leadtime} min",
+        title=f"Linear + {leadtime} min",
+        axis="off",
+        colorbar=False,
+    )
+
+    # Nowcast with salient blending into NWP
+    plt.subplot(n_leadtimes, 4, n * 4 + 3)
+    plot_precip_field(
+        precip_salient_blended[int(leadtime / timestep) - 1, :, :],
+        geodata=radar_metadata,
+        title=f"Salient + {leadtime} min",
         axis="off",
         colorbar=False,
     )
 
     # Raw NWP forecast
-    plt.subplot(n_leadtimes, 2, n * 2 + 2)
+    plt.subplot(n_leadtimes, 4, n * 4 + 4)
     plot_precip_field(
         precip_nwp[int(leadtime / timestep) - 1, :, :],
         geodata=nwp_metadata,
-        title=f"NWP +{leadtime} min",
+        title=f"NWP + {leadtime} min",
         axis="off",
         colorbar=False,
     )
 
 plt.tight_layout()
 plt.show()
+
+################################################################################
+# Note that the NaN values of the extrapolation forecast are replaced with NWP data
+# in the blended forecast, even before the blending starts.

pysteps/blending/interface.py

@@ -11,11 +11,14 @@
     get_method
 """
 
+from functools import partial
+
 from pysteps.blending import linear_blending
 from pysteps.blending import steps
 
 _blending_methods = dict()
 _blending_methods["linear_blending"] = linear_blending.forecast
+_blending_methods["salient_blending"] = partial(linear_blending.forecast, saliency=True)
 _blending_methods["steps"] = steps.forecast
 
 
@@ -32,6 +35,13 @@ def get_method(name):
     | linear_blending  | the linear blending of a nowcast method with other   |
     |                  | data (e.g. NWP data).                                |
     +------------------+------------------------------------------------------+
+    | salient_blending | the salient blending of a nowcast method with other  |
+    |                  | data (e.g. NWP data) described in :cite:`Hwang2015`. |
+    |                  | The blending is based on intensities and forecast    |
+    |                  | times. The blended product preserves pixel           |
+    |                  | intensities with time if they are strong enough based|
+    |                  | on their ranked salience.                            |
+    +------------------+------------------------------------------------------+
     | steps            | the STEPS stochastic nowcasting blending method      |
     |                  | described in :cite:`Seed2003`, :cite:`BPS2006` and   |
     |                  | :cite:`SPN2013`. The blending weights approach       |

pysteps/blending/linear_blending.py

@@ -9,9 +9,10 @@
 consists of the nowcast(s). In between the start and end time, the nowcast(s)
 weight decreases and NWP forecasts weight increases linearly from 1(0) to
 0(1). After the end time, the blended forecast entirely consists of the NWP
-forecasts.
+forecasts. The saliency-based blending method also takes into account the pixel
+intensities and preserves them if they are strong enough based on their ranked salience.
 
-Implementation of the linear blending between nowcast and NWP data.
+Implementation of the linear blending and saliency-based blending between nowcast and NWP data.
 
 .. autosummary::
     :toctree: ../generated/
@@ -22,6 +23,7 @@
 import numpy as np
 from pysteps import nowcasts
 from pysteps.utils import conversion
+from scipy.stats import rankdata
 
 
 def forecast(
@@ -36,10 +38,11 @@ def forecast(
     start_blending=120,
     end_blending=240,
     fill_nwp=True,
+    saliency=False,
     nowcast_kwargs=None,
 ):
 
-    """Generate a forecast by linearly blending nowcasts with NWP data
+    """Generate a forecast by linearly or saliency-based blending of nowcasts with NWP data
 
     Parameters
     ----------
@@ -81,12 +84,18 @@ def forecast(
     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.
+    saliency: bool, optional
+      Default value is False. If True, saliency will be used for blending. The blending
+      is based on intensities and forecast times as described in :cite:`Hwang2015`. The blended
+      product preserves pixel intensities with time if they are strong enough based on their ranked
+      salience.
     nowcast_kwargs: dict, optional
       Dictionary containing keyword arguments for the nowcast method.
 
+
     Returns
     -------
-    R_blended: ndarray
+    precip_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
@@ -166,45 +175,139 @@ def forecast(
         )
 
         # Initialise output
-        R_blended = np.zeros_like(precip_nowcast)
+        precip_blended = np.zeros_like(precip_nowcast)
 
         # Calculate the weights
         for i in range(timesteps):
             # Calculate what time we are at
             t = (i + 1) * timestep
 
+            if n_ens_members_max == 1:
+                ref_dim = 0
+            else:
+                ref_dim = 1
+
+            # apply blending
+            # compute the slice indices
+            slc_id = _get_slice(precip_blended.ndim, ref_dim, i)
+
             # 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:
+            if weight_nwp <= 0.0:
                 weight_nwp = 0.0
-            elif weight_nwp > 1.0:
-                weight_nwp = 1.0
+                precip_blended[slc_id] = precip_nowcast[slc_id]
 
-            # Calculate weight_nowcast
-            weight_nowcast = 1.0 - weight_nwp
+            elif weight_nwp >= 1.0:
+                weight_nwp = 1.0
+                precip_blended[slc_id] = precip_nwp[slc_id]
 
-            # 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, :, :]
-                )
+                # 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 saliency:
+                    ranked_salience = _get_ranked_salience(
+                        precip_nowcast[slc_id], precip_nwp[slc_id]
+                    )
+                    ws = _get_ws(weight_nowcast, ranked_salience)
+                    precip_blended[slc_id] = (
+                        ws * precip_nowcast[slc_id] + (1 - ws) * precip_nwp[slc_id]
+                    )
+
+                else:
+                    precip_blended[slc_id] = (
+                        weight_nwp * precip_nwp[slc_id]
+                        + weight_nowcast * precip_nowcast[slc_id]
+                    )
 
             # 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]
+                nan_indices = np.isnan(precip_blended)
+                precip_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
+        precip_blended = precip_nowcast
+
+    return precip_blended
+
+
+def _get_slice(n_dims, ref_dim, ref_id):
+    """source: https://stackoverflow.com/a/24399139/4222370"""
+    slc = [slice(None)] * n_dims
+    slc[ref_dim] = ref_id
+    return tuple(slc)
+
+
+def _get_ranked_salience(precip_nowcast, precip_nwp):
+    """Calculate ranked salience, which show how close the pixel is to the maximum intensity difference [r(x,y)=1]
+      or the minimum intensity difference [r(x,y)=0]
+
+    Parameters
+    ----------
+    precip_nowcast: array_like
+      Array of shape (m,n) containing the extrapolated precipitation field at a specified timestep
+    precip_nwp: array_like
+      Array of shape (m,n) containing the NWP fields at a specified timestep
+
+    Returns
+    -------
+    ranked_salience:
+      Array of shape (m,n) containing ranked salience
+    """
+
+    # calcutate normalized intensity
+    if np.max(precip_nowcast) == 0:
+        norm_nowcast = np.zeros_like(precip_nowcast)
+    else:
+        norm_nowcast = precip_nowcast / np.max(precip_nowcast)
+
+    if np.max(precip_nwp) == 0:
+        norm_nwp = np.zeros_like(precip_nwp)
+    else:
+        norm_nwp = precip_nwp / np.max(precip_nwp)
+
+    diff = norm_nowcast - norm_nwp
+
+    # Calculate ranked salience, based on dense ranking method, in which equally comparable values receive the same ranking number
+    ranked_salience = rankdata(diff, method="dense").reshape(diff.shape).astype("float")
+    ranked_salience /= ranked_salience.max()
+
+    return ranked_salience
+
 
-    return R_blended
+def _get_ws(weight, ranked_salience):
+    """Calculate salience weight based on linear weight and ranked salience as described in :cite:`Hwang2015`.
+      Cells with higher intensities result in larger weights.
+
+    Parameters
+    ----------
+    weight: int
+      Varying between 0 and 1
+    ranked_salience: array_like
+      Array of shape (m,n) containing ranked salience
+
+    Returns
+    -------
+    ws: array_like
+      Array of shape (m,n) containing salience weight, which preserves pixel intensties with time if they are strong
+      enough based on the ranked salience.
+    """
+
+    # Calculate salience weighte
+    ws = 0.5 * (
+        (weight * ranked_salience)
+        / (weight * ranked_salience + (1 - weight) * (1 - ranked_salience))
+        + (
+            np.sqrt(ranked_salience**2 + weight**2)
+            / (
+                np.sqrt(ranked_salience**2 + weight**2)
+                + np.sqrt((1 - ranked_salience) ** 2 + (1 - weight) ** 2)
+            )
+        )
+    )
+    return ws