Ide nowcasting (#221)

* Initial commit

* Docstring update

* Module name change

* Update module docstring

* Docstring additions

* Set interpolation order to 1

* Change interpolation order back to 0

* Add references to IDE-based nowcasting methods

* Implement utility methods

* Implement computation of convolution kernel

* Implement separate methods for isotropic and anisotropic kernels

* Implement optimization of convolution parameters

* Add the num_timesteps argument

* Implement computation of window weights

* Docstring revision

* Add missing imports

* Complete the input arguments and docstring

* New version of _optimize_convol_params

* New version of AR parameter optimization

* Implement iteration of the AR model

* Update the description of scikit-image

* Add LINDA reference

* Add LINDA to the nowcast interface

* First version for deterministic nowcasts

* Several fixes

* Add more printing and measurement of computation time

* Measure and return initialization and forecast loop times

* Fix estimation of ARI(2,1) parameters

* Remove unnecessary line

* Fix estimation of isotropic kernels

* Sort methods in alphabetical order

* Implement separate methods for initialization and computation of deterministic nowcasts

* Preparatory work for ensemble nowcasts

* Implement perturbation generator

* Add perturbations to deterministic nowcasts

* Implement ensemble nowcasts

* Implement seeding of the perturbation generator

* Add TODO comment

* Remove unimplemented option

* Arrange methods in alphabetical order

* Simplify convolution kernel and ACF computation

* Numerical integration adjustments

* Use exponential ACF instead of Gaussian

* Raise exception if non-integer timesteps are given

* Simplify parameter selection by using only one localization window radius for the deterministic model

* Add arguments and docstrings needed for velocity perturbations

* Minor refactoring and docstring revisions

* Add default values to docstring

* Add more printing

* Remove unnecessary arguments in _linda_init

* Implement advection field perturbation

* Refactor perturbation generators

* Allow input fields with odd dimensions

* Remove extra output argument

* Add default ARI model order to docstring

* Raise exception if no features are detected

* Add tests for LINDA nowcasts

* Make thresholds for perturbator estimation user-configurable

* Remove TODO comments

* Disable printing of warnings

* Add linda module

* Extend testing

* Limit number of realizations

* Upscale mm/h; log transform as an option

* Black formatting

* Adapt max crps for smaller ens size

* Specify in the docstring that only integer timesteps are accepted

* Remove unnecessary argument

* Use None instead of empty dict argument

* Disable test for fractional time steps

* Use f-strings instead of format

* Use lowercase variable names

* Simplify summation loop

* Use enumerate to simplify for loop

* Disable velocity perturbations by default

* Replace empty dictionary argument with None

* Simplify syntax of for loop

* Write all function documentation into docstrings

* Minor fixes

Comments, unused variables, and so on

* Implement measurement of coomputation time

* Use smaller domain for tests

* Black formatting

* Fix test

* Add note about the recommented number of features

* Remove unused arguments

* Test parallel computation with Dask

* Implement separate functions for initialization of the deterministic nowcast and the perturbation generator

* Docstring addition

* Add recommendation about using multiple parallel workers

* Use threads scheduler instead of multiprocessing

* Format with black

* Reduce the default number of ensemble members to 10

* Add option to use multiprocessing

* Add max_num_features argument

* First version of LINDA examples

* Change variable name

* Fix incorrect dictionary names

* Disable warnings

* Remove blank space before colons in docstrings

* Plot S-PROG nowcast for comparison

* Add discussion about the LINDA/S-PROG comparison

* Add comparison with STEPS

* Adjust figure size

* Remove unused import

* Adjust title text

* Change file name

* Remove ticks

Co-authored-by: Seppo Pulkkinen <[email protected]>
Co-authored-by: ned <[email protected]>

Author: 3 people

doc/source/pysteps_reference/nowcasts.rst

@@ -8,6 +8,7 @@ Implementation of deterministic and ensemble nowcasting methods.
 .. automodule:: pysteps.nowcasts.interface
 .. automodule:: pysteps.nowcasts.anvil
 .. automodule:: pysteps.nowcasts.extrapolation
+.. automodule:: pysteps.nowcasts.linda
 .. automodule:: pysteps.nowcasts.lagrangian_probability
 .. automodule:: pysteps.nowcasts.sprog
 .. automodule:: pysteps.nowcasts.sseps

doc/source/references.bib

@@ -73,6 +73,16 @@ @ARTICLE{FNPC2020
   DOI = "10.3390/atmos11030267"
 }
 
+@ARTICLE{FW2005,
+  AUTHOR = "N. I. Fox and C. K. Wikle",
+  TITLE = "A Bayesian Quantitative Precipitation Nowcast Scheme",
+  JOURNAL = "Weather and Forecasting",
+  VOLUME = 20,
+  NUMBER = 3,
+  PAGES = "264--275",
+  YEAR = 2005
+}
+
 @ARTICLE{GZ2002,
   AUTHOR = "U. Germann and I. Zawadzki",
   TITLE = "Scale-Dependence of the Predictability of Precipitation from Continental Radar Images. {P}art {I}: Description of the Methodology",
@@ -157,6 +167,14 @@ @ARTICLE{PCLH2020
   YEAR = 2020
 }
 
+@ARTICLE{PCN2021,
+  AUTHOR = "S. Pulkkinen and V. Chandrasekar and T. Niemi",
+  TITLE = "Lagrangian Integro-Difference Equation Model for Precipitation Nowcasting",
+  JOURNAL = "Journal of Atmospheric and Oceanic Technology",
+  NOTE = "submitted",
+  YEAR = 2021
+}
+
 @INCOLLECTION{PGPO1994,
   AUTHOR = "M. Proesmans and L. van Gool and E. Pauwels and A. Oosterlinck",
   TITLE = "Determination of optical flow and its discontinuities using non-linear diffusion",
@@ -235,6 +253,16 @@ @ARTICLE{SPN2013
   DOI = "10.1002/wrcr.20536"
 }
 
+@ARTICLE{XWF2005,
+  AUTHOR = "K. Xu and C. K Wikle and N. I. Fox",
+  TITLE = "A Kernel-Based Spatio-Temporal Dynamical Model for Nowcasting Weather Radar Reflectivities",
+  JOURNAL = "Journal of the American Statistical Association",
+  VOLUME = 100,
+  NUMBER = 472,
+  PAGES = "1133--1144",
+  YEAR  = 2005
+}
+
 @ARTICLE{ZR2009,
   AUTHOR = "P. Zacharov and D. Rezacova",
   TITLE = "Using the fractions skill score to assess the relationship between an ensemble {QPF} spread and skill",

doc/source/user_guide/install_pysteps.rst

@@ -36,7 +36,7 @@ Other optional dependencies include:
 * `pywavelets <https://pywavelets.readthedocs.io/en/latest/>`_
   (for intensity-scale verification)
 * `pandas <https://pandas.pydata.org/>`_ and
-  `scikit-image <https://scikit-image.org/>`_ (for more advanced feature tracking)
+  `scikit-image <https://scikit-image.org/>`_ (for the DATing and LINDA nowcast methods)
 
 **Important**: If you only want to use pysteps, you can continue reading below.
 But, if you want to contribute to pysteps or edit the package, you need to install

examples/linda_nowcasts.py

@@ -0,0 +1,185 @@
+#!/bin/env python
+"""
+LINDA nowcasts
+==============
+
+This example shows how to compute and plot a deterministic and ensemble LINDA
+nowcasts using Swiss radar data.
+
+"""
+
+from datetime import datetime
+import warnings
+
+warnings.simplefilter("ignore")
+
+import matplotlib.pyplot as plt
+
+from pysteps import io, rcparams
+from pysteps.motion.lucaskanade import dense_lucaskanade
+from pysteps.nowcasts import linda, sprog, steps
+from pysteps.utils import conversion, dimension, transformation
+from pysteps.visualization import plot_precip_field
+
+###############################################################################
+# Read the input rain rate fields
+# -------------------------------
+
+date = datetime.strptime("201701311200", "%Y%m%d%H%M")
+data_source = "mch"
+
+# Read the data source information from rcparams
+datasource_params = rcparams.data_sources[data_source]
+
+# Find the radar files in the archive
+fns = io.find_by_date(
+    date,
+    datasource_params["root_path"],
+    datasource_params["path_fmt"],
+    datasource_params["fn_pattern"],
+    datasource_params["fn_ext"],
+    datasource_params["timestep"],
+    num_prev_files=2,
+)
+
+# Read the data from the archive
+importer = io.get_method(datasource_params["importer"], "importer")
+reflectivity, _, metadata = io.read_timeseries(
+    fns, importer, **datasource_params["importer_kwargs"]
+)
+
+# Convert reflectivity to rain rate
+rainrate, metadata = conversion.to_rainrate(reflectivity, metadata)
+
+# Upscale data to 2 km to reduce computation time
+rainrate, metadata = dimension.aggregate_fields_space(rainrate, metadata, 2000)
+
+# Plot the most recent rain rate field
+plt.figure()
+plot_precip_field(rainrate[-1, :, :])
+plt.show()
+
+###############################################################################
+# Estimate the advection field
+# ----------------------------
+
+# The advection field is estimated using the Lucas-Kanade optical flow
+advection = dense_lucaskanade(rainrate, verbose=True)
+
+###############################################################################
+# Deterministic nowcast
+# ---------------------
+
+# Compute 30-minute LINDA nowcast with 8 parallel workers
+# Restrict the number of features to 15 to reduce computation time
+nowcast_linda = linda.forecast(
+    rainrate,
+    advection,
+    6,
+    max_num_features=15,
+    add_perturbations=False,
+    num_workers=8,
+    measure_time=True,
+)[0]
+
+# Compute S-PROG nowcast for comparison
+rainrate_db, _ = transformation.dB_transform(
+    rainrate, metadata, threshold=0.1, zerovalue=-15.0
+)
+nowcast_sprog = sprog.forecast(
+    rainrate_db[-3:, :, :],
+    advection,
+    6,
+    n_cascade_levels=6,
+    R_thr=-10.0,
+)
+
+# Convert reflectivity nowcast to rain rate
+nowcast_sprog = transformation.dB_transform(
+    nowcast_sprog, threshold=-10.0, inverse=True
+)[0]
+
+# Plot the nowcasts
+fig = plt.figure(figsize=(9, 4))
+ax = fig.add_subplot(1, 2, 1)
+plot_precip_field(
+    nowcast_linda[-1, :, :],
+    title="LINDA (+ 30 min)",
+)
+
+ax = fig.add_subplot(1, 2, 2)
+plot_precip_field(
+    nowcast_sprog[-1, :, :],
+    title="S-PROG (+ 30 min)",
+)
+
+plt.show()
+
+###############################################################################
+# The above figure shows that the filtering scheme implemented in LINDA preserves
+# small-scale and band-shaped features better than S-PROG. This is because the
+# former uses a localized elliptical convolution kernel instead of the
+# cascade-based autoregressive process, where the parameters are estimated over
+# the whole domain.
+
+###############################################################################
+# Probabilistic nowcast
+# ---------------------
+
+# Compute 30-minute LINDA nowcast ensemble with 40 members and 8 parallel workers
+nowcast_linda = linda.forecast(
+    rainrate,
+    advection,
+    6,
+    max_num_features=15,
+    add_perturbations=True,
+    num_ens_members=40,
+    num_workers=8,
+    measure_time=True,
+)[0]
+
+# Compute 40-member STEPS nowcast for comparison
+nowcast_steps = steps.forecast(
+    rainrate_db[-3:, :, :],
+    advection,
+    6,
+    40,
+    n_cascade_levels=6,
+    R_thr=-10.0,
+    mask_method="incremental",
+    kmperpixel=2.0,
+    timestep=datasource_params["timestep"],
+    vel_pert_method=None,
+)
+
+# Convert reflectivity nowcast to rain rate
+nowcast_steps = transformation.dB_transform(
+    nowcast_steps, threshold=-10.0, inverse=True
+)[0]
+
+# Plot two ensemble members of both nowcasts
+fig = plt.figure()
+for i in range(2):
+    ax = fig.add_subplot(2, 2, i + 1)
+    ax = plot_precip_field(
+        nowcast_linda[i, -1, :, :], geodata=metadata, colorbar=False, axis="off"
+    )
+    ax.set_title(f"LINDA Member {i+1}")
+
+for i in range(2):
+    ax = fig.add_subplot(2, 2, 3 + i)
+    ax = plot_precip_field(
+        nowcast_steps[i, -1, :, :], geodata=metadata, colorbar=False, axis="off"
+    )
+    ax.set_title(f"STEPS Member {i+1}")
+
+###############################################################################
+# The above figure shows the main difference between LINDA and STEPS. In
+# addition to the convolution kernel, another improvement in LINDA is a
+# localized perturbation generator using the short-space Fourier transform
+# (SSFT) and a spatially variable marginal distribution. As a result, the
+# LINDA ensemble members preserve the anisotropic and small-scale structures
+# considerably better than STEPS.
+
+plt.tight_layout()
+plt.show()

pysteps/nowcasts/interface.py

@@ -32,16 +32,17 @@
 """
 
 from pysteps.extrapolation.interface import eulerian_persistence
-from pysteps.nowcasts import anvil, sprog, steps, sseps, extrapolation
+from pysteps.nowcasts import anvil, extrapolation, linda, sprog, steps, sseps
 from pysteps.nowcasts import lagrangian_probability
 
 _nowcast_methods = dict()
 _nowcast_methods["anvil"] = anvil.forecast
 _nowcast_methods["eulerian"] = eulerian_persistence
 _nowcast_methods["extrapolation"] = extrapolation.forecast
 _nowcast_methods["lagrangian"] = extrapolation.forecast
-_nowcast_methods["probability"] = lagrangian_probability.forecast
 _nowcast_methods["lagrangian_probability"] = lagrangian_probability.forecast
+_nowcast_methods["linda"] = linda.forecast
+_nowcast_methods["probability"] = lagrangian_probability.forecast
 _nowcast_methods["sprog"] = sprog.forecast
 _nowcast_methods["sseps"] = sseps.forecast
 _nowcast_methods["steps"] = steps.forecast
@@ -68,8 +69,10 @@ def get_method(name):
     |  lagrangian or  | this approach extrapolates the last observation       |
     |  extrapolation  | using the motion field (Lagrangian persistence)       |
     +-----------------+-------------------------------------------------------+
+    |  linda          | the LINDA method developed in :cite:`PCN2021`         |
+    +-----------------+-------------------------------------------------------+
     |  lagrangian_\   | this approach computes local lagrangian probability   |
-    |  probability    | forecasts of threshold exceedences.                   |
+    |  probability    | forecasts of threshold exceedences                    |
     +-----------------+-------------------------------------------------------+
     |  sprog          | the S-PROG method described in :cite:`Seed2003`       |
     +-----------------+-------------------------------------------------------+
@@ -78,7 +81,7 @@ def get_method(name):
     |                 |                                                       |
     +-----------------+-------------------------------------------------------+
     |  sseps          | short-space ensemble prediction system (SSEPS).       |
-    |                 | Essentially, this is a localization of STEPS.         |
+    |                 | Essentially, this is a localization of STEPS          |
     +-----------------+-------------------------------------------------------+
     """
     if isinstance(name, str):