Steps blending v1

.github/workflows/test_pysteps.yml

@@ -3,9 +3,13 @@ name: Test pysteps
 on:
   # Triggers the workflow on push or pull request events to the master branch
   push:
-    branches: [ master ]
+    branches:
+      - master
+      - pysteps-v2
   pull_request:
-    branches: [ master ]
+      branches:
+        - master
+        - pysteps-v2
 
 jobs:
   unit_tests:

.pre-commit-config.yaml

@@ -1,6 +1,6 @@
 repos:
 -   repo: https://github.com/psf/black
-    rev: 21.6b0
+    rev: 21.7b0
     hooks:
     - id: black
       language_version: python3

ci/ci_test_env.yml

@@ -26,6 +26,7 @@ dependencies:
   - PyWavelets
   - pandas
   - scikit-image
+  - rasterio
   # Test dependencies
   - pytest
   - pytest-cov

doc/requirements.txt

@@ -1,8 +1,8 @@
-# Additional requeriments related to the documentation build only
+# Additional requirements related to the documentation build only
 sphinx
 sphinxcontrib.bibtex
 sphinx-book-theme
 sphinx_gallery
 scikit-image
 pandas
-
+git+https://github.com/pySTEPS/pysteps-nwp-importers.git@main#egg=pysteps_nwp_importers

doc/source/pysteps_reference/blending.rst

@@ -0,0 +1,12 @@
+================
+pysteps.blending
+================
+
+Implementation of blending methods for blending (ensemble) nowcasts with Numerical Weather Prediction (NWP) models.
+
+.. automodule:: pysteps.blending.interface
+.. automodule:: pysteps.blending.clim
+.. automodule:: pysteps.blending.linear_blending
+.. automodule:: pysteps.blending.skill_scores
+.. automodule:: pysteps.blending.steps
+.. automodule:: pysteps.blending.utils

doc/source/pysteps_reference/index.rst

@@ -14,6 +14,7 @@ available in pysteps.
     :caption: API Reference
 
     pysteps
+    blending
     cascade
     decorators
     extrapolation

doc/source/pysteps_reference/utils.rst

@@ -16,3 +16,4 @@ Implementation of miscellaneous utility functions.
 .. automodule:: pysteps.utils.spectral
 .. automodule:: pysteps.utils.tapering
 .. automodule:: pysteps.utils.transformation
+.. automodule:: pysteps.utils.reprojection

doc/source/references.bib

@@ -1,4 +1,14 @@
 
+@TECHREPORT{BPS2004,
+  AUTHOR = "N. E. Bowler and C. E. Pierce and A. W. Seed",
+  TITLE = "{STEPS}: A probabilistic precipitation forecasting scheme which merges an extrapolation nowcast with downscaled {NWP}",
+  INSTITUTION = "UK Met Office",
+  TYPE = "Forecasting Research Technical Report",
+  NUMBER = 433,
+  ADDRESS = "Wallingford, United Kingdom",
+  YEAR = 2004,
+}
+
 @ARTICLE{BPS2006,
   AUTHOR = "N. E. Bowler and C. E. Pierce and A. W. Seed",
   TITLE = "{STEPS}: A probabilistic precipitation forecasting scheme which merges an extrapolation nowcast with downscaled {NWP}",

doc/source/user_guide/install_pysteps.rst

@@ -36,7 +36,9 @@ 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 the DATing and LINDA nowcast methods)
+  `scikit-image <https://scikit-image.org/>`_ (for advanced feature detection methods)
+* `rasterio <https://rasterio.readthedocs.io/en/latest/>`_ (for the reprojection module)
+
 
 **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

environment_dev.yml

@@ -30,3 +30,4 @@ dependencies:
   - cartopy>=0.18
   - scikit-image
   - pandas
+  - rasterio

examples/anvil_nowcast.py

@@ -118,9 +118,7 @@
     date, root_path, path_fmt, fn_pattern, fn_ext, timestep=5, num_next_files=3
 )
 
-refobs_field, quality, metadata = io.read_timeseries(
-    filenames, importer, **importer_kwargs
-)
+refobs_field, _, metadata = io.read_timeseries(filenames, importer, **importer_kwargs)
 
 refobs_field, metadata = utils.to_rainrate(refobs_field[-1], metadata)
 refobs_field[refobs_field < 0.5] = 0.0

examples/blended_forecast.py

@@ -0,0 +1,230 @@
+# -*- coding: utf-8 -*-
+"""
+Blended forecast
+====================
+
+This tutorial shows how to construct a blended forecast from an ensemble nowcast
+using the STEPS approach and a Numerical Weather Prediction (NWP) rainfall
+forecast. The used datasets are from the Bureau of Meteorology, Australia.
+
+"""
+
+import os
+from datetime import datetime
+
+import numpy as np
+from matplotlib import pyplot as plt
+
+import pysteps
+from pysteps import io, rcparams, blending
+from pysteps.visualization import plot_precip_field
+
+
+################################################################################
+# Read the radar images and the NWP forecast
+# ------------------------------------------
+#
+# First, we import a sequence of 3 images of 10-minute radar composites
+# and the corresponding NWP rainfall forecast that was available at that time.
+#
+# You need the pysteps-data archive downloaded and the pystepsrc file
+# configured with the data_source paths pointing to data folders.
+# Additionally, the pysteps-nwp-importers plugin needs to be installed, see
+# https://github.com/pySTEPS/pysteps-nwp-importers.
+
+# Selected case
+date_radar = datetime.strptime("202010310400", "%Y%m%d%H%M")
+# The last NWP forecast was issued at 00:00
+date_nwp = datetime.strptime("202010310000", "%Y%m%d%H%M")
+radar_data_source = rcparams.data_sources["bom"]
+nwp_data_source = rcparams.data_sources["bom_nwp"]
+
+###############################################################################
+# Load the data from the archive
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+root_path = radar_data_source["root_path"]
+path_fmt = "prcp-c10/66/%Y/%m/%d"
+fn_pattern = "66_%Y%m%d_%H%M00.prcp-c10"
+fn_ext = radar_data_source["fn_ext"]
+importer_name = radar_data_source["importer"]
+importer_kwargs = radar_data_source["importer_kwargs"]
+timestep = 10.0
+
+# Find the radar files in the archive
+fns = io.find_by_date(
+    date_radar, root_path, path_fmt, fn_pattern, fn_ext, timestep, num_prev_files=2
+)
+
+# Read the radar composites
+importer = io.get_method(importer_name, "importer")
+radar_precip, _, radar_metadata = io.read_timeseries(fns, importer, **importer_kwargs)
+
+# Import the NWP data
+filename = os.path.join(
+    nwp_data_source["root_path"],
+    datetime.strftime(date_nwp, nwp_data_source["path_fmt"]),
+    datetime.strftime(date_nwp, nwp_data_source["fn_pattern"])
+    + "."
+    + nwp_data_source["fn_ext"],
+)
+
+nwp_importer = io.get_method("bom_nwp", "importer")
+nwp_precip, _, nwp_metadata = nwp_importer(filename)
+
+# Only keep the NWP forecasts from the last radar observation time (2020-10-31 04:00)
+# onwards
+
+nwp_precip = nwp_precip[24:43, :, :]
+
+
+################################################################################
+# Pre-processing steps
+# --------------------
+
+# Make sure the units are in mm/h
+converter = pysteps.utils.get_method("mm/h")
+radar_precip, radar_metadata = converter(radar_precip, radar_metadata)
+nwp_precip, nwp_metadata = converter(nwp_precip, nwp_metadata)
+
+# Threshold the data
+radar_precip[radar_precip < 0.1] = 0.0
+nwp_precip[nwp_precip < 0.1] = 0.0
+
+# Plot the radar rainfall field and the first time step of the NWP forecast.
+date_str = datetime.strftime(date_radar, "%Y-%m-%d %H:%M")
+plt.figure(figsize=(10, 5))
+plt.subplot(121)
+plot_precip_field(
+    radar_precip[-1, :, :],
+    geodata=radar_metadata,
+    title=f"Radar observation at {date_str}",
+)
+plt.subplot(122)
+plot_precip_field(
+    nwp_precip[0, :, :], geodata=nwp_metadata, title=f"NWP forecast at {date_str}"
+)
+plt.tight_layout()
+plt.show()
+
+# transform the data to dB
+transformer = pysteps.utils.get_method("dB")
+radar_precip, radar_metadata = transformer(radar_precip, radar_metadata, threshold=0.1)
+nwp_precip, nwp_metadata = transformer(nwp_precip, nwp_metadata, threshold=0.1)
+
+# r_nwp has to be four dimentional (n_models, time, y, x).
+# If we only use one model:
+if nwp_precip.ndim == 3:
+    nwp_precip = nwp_precip[None, :]
+
+###############################################################################
+# For the initial time step (t=0), the NWP rainfall forecast is not that different
+# from the observed radar rainfall, but it misses some of the locations and
+# shapes of the observed rainfall fields. Therefore, the NWP rainfall forecast will
+# initially get a low weight in the blending process.
+#
+# Determine the velocity fields
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+oflow_method = pysteps.motion.get_method("lucaskanade")
+
+# First for the radar images
+velocity_radar = oflow_method(radar_precip)
+
+# Then for the NWP forecast
+velocity_nwp = []
+# Loop through the models
+for n_model in range(nwp_precip.shape[0]):
+    # Loop through the timesteps. We need two images to construct a motion
+    # field, so we can start from timestep 1. Timestep 0 will be the same
+    # as timestep 1.
+    _v_nwp_ = []
+    for t in range(1, nwp_precip.shape[1]):
+        v_nwp_ = oflow_method(nwp_precip[n_model, t - 1 : t + 1, :])
+        _v_nwp_.append(v_nwp_)
+        v_nwp_ = None
+    # Add the velocity field at time step 1 to time step 0.
+    _v_nwp_ = np.insert(_v_nwp_, 0, _v_nwp_[0], axis=0)
+    velocity_nwp.append(_v_nwp_)
+velocity_nwp = np.stack(velocity_nwp)
+
+
+################################################################################
+# The blended forecast
+# --------------------
+
+precip_forecast = blending.steps.forecast(
+    precip=radar_precip,
+    precip_models=nwp_precip,
+    velocity=velocity_radar,
+    velocity_models=velocity_nwp,
+    timesteps=18,
+    timestep=timestep,
+    issuetime=date_radar,
+    n_ens_members=1,
+    precip_thr=radar_metadata["threshold"],
+    kmperpixel=radar_metadata["xpixelsize"] / 1000.0,
+    noise_stddev_adj="auto",
+    vel_pert_method=None,
+)
+
+# Transform the data back into mm/h
+precip_forecast, _ = converter(precip_forecast, radar_metadata)
+radar_precip, _ = converter(radar_precip, radar_metadata)
+nwp_precip, _ = converter(nwp_precip, nwp_metadata)
+
+
+################################################################################
+# Visualize the output
+# ~~~~~~~~~~~~~~~~~~~~
+#
+# The NWP rainfall forecast has a lower weight than the radar-based extrapolation
+# forecast at the issue time of the forecast (+0 min). Therefore, the first time
+# steps consist mostly of the extrapolation.
+# However, near the end of the forecast (+180 min), the NWP share in the blended
+# forecast has become more important and the forecast starts to resemble the
+# NWP forecast more.
+
+fig = plt.figure(figsize=(5, 12))
+
+leadtimes_min = [30, 60, 90, 120, 150, 180]
+n_leadtimes = len(leadtimes_min)
+for n, leadtime in enumerate(leadtimes_min):
+
+    # Nowcast with blending into NWP
+    plt.subplot(n_leadtimes, 2, n * 2 + 1)
+    plot_precip_field(
+        precip_forecast[0, int(leadtime / timestep) - 1, :, :],
+        geodata=radar_metadata,
+        title=f"Nowcast +{leadtime} min",
+        axis="off",
+        colorbar=False,
+    )
+
+    # Raw NWP forecast
+    plt.subplot(n_leadtimes, 2, n * 2 + 2)
+    plot_precip_field(
+        nwp_precip[0, int(leadtime / timestep) - 1, :, :],
+        geodata=nwp_metadata,
+        title=f"NWP +{leadtime} min",
+        axis="off",
+        colorbar=False,
+    )
+
+
+################################################################################
+# References
+# ~~~~~~~~~~
+#
+# Bowler, N. E., and C. E. Pierce, and A. W. Seed. 2004. "STEPS: A probabilistic
+# precipitation forecasting scheme which merges an extrapolation nowcast with
+# downscaled NWP." Forecasting Research Technical Report No. 433. Wallingford, UK.
+#
+# Bowler, N. E., and C. E. Pierce, and A. W. Seed. 2006. "STEPS: A probabilistic
+# precipitation forecasting scheme which merges an extrapolation nowcast with
+# downscaled NWP." Quarterly Journal of the Royal Meteorological Society 132(16):
+# 2127-2155. https://doi.org/10.1256/qj.04.100
+#
+# Seed, A. W., and C. E. Pierce, and K. Norman. 2013. "Formulation and evaluation
+# of a scale decomposition-based stochastic precipitation nowcast scheme." Water
+# Resources Research 49(10): 6624-664. https://doi.org/10.1002/wrcr.20536