Implement STEPS blending (#233)

* First basic functions to implement STEPS blending * Add compute of blend means,sigmas and recompose * pysteps.io with xarray (#219) * Add xarray dependency * MCH importer returns an xarray Dataset * Remove plot lines * Remove import * Adapt readers to xarray format * Rewrite as more general decorator * Add missing metadata * Adapt io tests * Mrms bounding box (#222) * Fix bounding box coordinates * Add missing metadata * Import xarray DataArray Ignore quality field * Black * Do not hardcode metadata * Address review comments by ruben * Add a legacy option to the io functions A "legacy" options is added to revert back the importers and readers behavior to version 1. This is a temporary solution to allow the examples, and other functions, to run as usual (v1.*). Hopefully, this is will allow a easier transition into version 2 during the development process and will allow testing functions that were not updated to v2. * Fix compatibility problems with tests Many of the tests were updated to use the legacy data structures (v1). The tests that still contains issues, were tagged with a TODO message and they are skipped. This will allow incremental changes to be tested in the new v2. IMPORTANT: once the v2 branch is stable, we may remove the legacy compatibility from the code base and the tests. * Update dependencies * Ignore plugins test Co-authored-by: Andres Perez Hortal <16256571+aperezhortal@users.noreply.github.com> * Add blend_optical_flow * changes to steps blending procedure - weights according to adjusted BPS2006 method * changes to blending procedures - adjust weights from original BPS2006 method * Determine spatial correlation of NWP model forecast * First attempt to make correlations and thus weights lead time dependent (in progress..) * Change back to original BPS2006 blending formulation and add regression of skill values to climatological values for weights determination * Reformat code with Black * Skill score script imports climatological correlation-values from file now * Small changes to skill score script * Add skill score tests and an interface * Add skill score tests and an interface * Small change to docstring * Bom import xarray (#228) * Add import_bom_rf3 using xarray * Add tests to xarray version * Fix mrms importer tests * Pass **kwargs to internal functions * Add nwp_importers to read bom nwp sample data * Add bom nwp data to source file * Add tests for bom_nwp reader * Fix pystepsrc Co-authored-by: Andres Perez Hortal <16256571+aperezhortal@users.noreply.github.com> * Functions to store and compute climatological weights (#231) * Implement the functions get_default_weights, save_weights, calc_clim_weights. These functions are used to evolve the weights in the scale- and skill-dependent blending with NWP in the STEPS blending algorithm. The current weights, based on the correlations per cascade level, are regressed towards these climatological weights in the course of the forecast. These functions save the current and compute the climatological weights (a running mean of the weights of the past n days, where typically n=30). First daily averages are stored and these are then averaged over the running window of n days. * Add tests for pysteps climatological weight io and calculations. * Add path_workdir to outputs section in pystepsrc file and use it as a default path to store/retrieve blending weights. * Minor changes to docstrings, changes to skill scores and testing scripts * Completed documentation for blending clim module, cleanup. Co-authored-by: RubenImhoff <r.o.imhoff@live.nl> * Main blending module, first steps * Add simple tests * Minor changes to tester: velocity now based on rainfall field of NWP * Add utilities to decompose, store and load NWP cascades for use in blending (#232) * First version of NWP decomposition * Added saving to netCDF * Completed functions for saving and loading decomposed NWP data * Added example files for the decomposed NWP functions * Added compatibility with numpy datetime * Use default output path_workdir for tmp files in blending/utils.py. * Update documentation of NWP decomposition functions in utils.py Co-authored-by: Wout Dewettinck <wout.dewettinck@ugent.be> Co-authored-by: wdewettin <87696913+wdewettin@users.noreply.github.com> * Add importer for RMI NWP data (#234) Add importer for netcdf NWP data from RMI using xarrays. * Add test for RMI NWP data importer. * Add entry for RMI NWP data in pystepsrc. * Run black on everything: fix formatting. * Add KNMI Harmonie NWP netcdf importer and tests (#235) * Changes to v_models to make use of multiple timesteps. Changes in the velocity field over time in the NWP forecast will be taken into account now. * Fixes for KNMI importer: Add forgotten @postprocess_import() Don't call dropna on NWP data. * Avoid shadowing of pysteps.blending.utils by pysteps.utils * First attempt for probability matching and masking utility; part 1 * Changes to prob matching and masking methods; part 2 * Prob matching and masking changes; part 3. Ready for testing with real data from here on * Remove unnecessary print statements * Cleanup imports * More cleanup * Update docstrings * RMI importer for gallery example (will follow) * Reprojection functionality (#236) * Added Lesley's reprojection module to this branch * Added compatibility for three-dimensional xarrays * Add commentary to reprojection util * Changes to make reprojection of KNMI data possible * Changes after Daniele's review * Add dependencies * Changes to importers, see issue #215 * Add tests * Fix some issues * documentation * Fixes for tests * Set requirements again * Some fixes * Changes to nwp_importers after Carlos' response * Remove wrong example script * Remove rasterio dependencies from lists * First try to prevent testing error * Changes Daniele and fix knmi nwp importer * Add rasterio to tox.ini * Aesthetics * rasterio import test * Add rasterio to the test dependencies * Reset try-except functionality for rasterio import * Fix for failing test on windows python 3.6 * add importerskip rasterio Co-authored-by: Wout Dewettinck <wout.dewettinck@ugent.be> * Fixes in nwp importers * Revert "Merge branch 'steps_blending' into pysteps-v2" (#239) This reverts commit 2c639f8, reversing changes made to bccb8fc. * Merge latest version pysteps-v2 into steps_blending branch (#237) * Update docstrings * More cleanup * Cleanup imports * Cleanup imports * More cleanup * Update docstrings * Update references Mention the work of Ravuri et al (2021, Nature) as an example of work using cGANs to generate ensembles * Clean up page * Reprojection functionality (#236) * Added Lesley's reprojection module to this branch * Added compatibility for three-dimensional xarrays * Add commentary to reprojection util * Changes to make reprojection of KNMI data possible * Changes after Daniele's review * Add dependencies * Changes to importers, see issue #215 * Add tests * Fix some issues * documentation * Fixes for tests * Set requirements again * Some fixes * Changes to nwp_importers after Carlos' response * Remove wrong example script * Remove rasterio dependencies from lists * First try to prevent testing error * Changes Daniele and fix knmi nwp importer * Add rasterio to tox.ini * Aesthetics * rasterio import test * Add rasterio to the test dependencies * Reset try-except functionality for rasterio import * Fix for failing test on windows python 3.6 * add importerskip rasterio Co-authored-by: Wout Dewettinck <wout.dewettinck@ugent.be> * Revert "Merge branch 'steps_blending' into pysteps-v2" (#239) This reverts commit 2c639f8, reversing changes made to bccb8fc. Co-authored-by: ned <daniele.nerini@meteoswiss.ch> Co-authored-by: dnerini <daniele.nerini@gmail.com> Co-authored-by: Wout Dewettinck <wout.dewettinck@ugent.be> * NWP skill calculation only within radar domain * Update docs * Add example for gallery examples * Fix docstrings example * Remove additional normalization step * Fixes for the tests * update docs * changes to post-processing rainfall field and docstrings * Update contributing guidelines (#241) - Improve grammar. - Make the guide more concise. Remove unused/unnecessary rules. - Indicate more clearly which parts of the guidelines are inspired by other projects (before they were only mentioned at the end). - Change "Travis-CI" references by "GitHub Actions". * Advect noise cascade * Allow for moving domain mask of extrapolation component * minor fixes * Linear blending (#229) * Implemented linear blending function * Added example file and test * Added compatibility for NWP ensembles The PR is ready to go. Making the code xarray ready will be done in a separate PR. Co-authored-by: RubenImhoff <r.o.imhoff@live.nl> * weights calculation adjustment outside radar domain if only one model present * allow for mirroring of advected noise cascade * implementation of weights following Seed et al. (2013) * Allow for decomposed NWP precip and NWP velocity fields: part 2 * Store decomposed fields with compression * changes after first review Daniele * Remove unnecessary print statement * fixes to blending utils and implementation of blending utils tests * remove unnecessary lines * Fix one time step shift of extrapolation skill prior to blending * minor changes to blending climatology, blending weights and remove path_workdir from pystepsrc * Make NWP forecast decomposition prior to blending function optional * Use pathlib * Extract methods * Minor changes to docstrings * Access climatological skill file for multiple NWP model and date string changes to prevent errors in blending.utils Co-authored-by: Carlos Velasco <carlos.velasco@bom.gov.au> Co-authored-by: ned <daniele.nerini@meteoswiss.ch> Co-authored-by: Andres Perez Hortal <16256571+aperezhortal@users.noreply.github.com> Co-authored-by: Ruben Imhoff <Ruben.Imhoff@deltares.nl> Co-authored-by: Carlos Velasco <cvelascof@gmail.com> Co-authored-by: Lesley De Cruz <lesley.decruz+git@gmail.com> Co-authored-by: Wout Dewettinck <wout.dewettinck@ugent.be> Co-authored-by: wdewettin <87696913+wdewettin@users.noreply.github.com> Co-authored-by: Lesley De Cruz <lesley.decruz@meteo.be> Co-authored-by: dnerini <daniele.nerini@gmail.com>
pySTEPS · Jan 14, 2022 · c71635a · c71635a
1 parent 82a9d18
commit c71635a
Show file tree

Hide file tree

Showing 27 changed files with 5,073 additions and 56 deletions.
diff --git a/doc/source/pysteps_reference/blending.rst b/doc/source/pysteps_reference/blending.rst
@@ -0,0 +1,14 @@
+================
+pysteps.blending
+================
+
+Implementation of blending methods for blending (ensemble) nowcasts with Numerical Weather Prediction (NWP) models.
+
+.. automodule:: pysteps.blending.clim
+.. automodule:: pysteps.blending.interface
+.. automodule:: pysteps.blending.skill_scores
+.. automodule:: pysteps.blending.utils
+.. automodule:: pysteps.blending.linear_blending
+.. automodule:: pysteps.blending.steps
+
+
diff --git a/doc/source/references.bib b/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}",

diff --git a/examples/blended_forecast.py b/examples/blended_forecast.py
@@ -0,0 +1,284 @@
+# -*- 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 Royal Meteorological Insitute of Belgium.  
+
+"""
+
+from matplotlib import pyplot as plt
+import numpy as np
+import os
+from datetime import datetime
+
+import pysteps
+from pysteps import io, rcparams, cascade, 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.
+
+# 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")
+r_radar = io.read_timeseries(fns, importer, legacy=False, **importer_kwargs)
+radar_data_xr = r_radar[-1, :, :]
+
+# Get the metadata
+radar_metadata = radar_data_xr.x.attrs.copy()
+radar_metadata.update(**radar_data_xr.y.attrs)
+radar_metadata.update(**radar_data_xr.attrs)
+
+# 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_data_xr = io.import_bom_nwp_xr(filename)
+nwp_metadata = nwp_data_xr.x.attrs.copy()
+nwp_metadata.update(**nwp_data_xr.y.attrs)
+nwp_metadata.update(**nwp_data_xr.attrs)
+
+# Only keep the NWP forecasts from the last radar observation time (2020-10-31 04:00)
+# onwards
+r_nwp = nwp_data_xr.sel(
+    t=slice(np.datetime64("2020-10-31T04:00"), np.datetime64("2020-10-31T07:00"))
+)
+
+
+################################################################################
+# Pre-processing steps
+# --------------------
+
+# Make sure the units are in mm/h
+converter = pysteps.utils.get_method("mm/h")
+r_radar, radar_metadata = converter(r_radar, radar_metadata)
+r_nwp, nwp_metadata = converter(r_nwp, nwp_metadata)
+
+# Threshold the data
+r_radar.data[r_radar.data < 0.1] = 0.0
+r_nwp.data[r_nwp.data < 0.1] = 0.0
+
+# Plot the radar rainfall field and the first time step of the NWP forecast.
+# 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.
+date_str = datetime.strftime(date_radar, "%Y-%m-%d %H:%M")
+plt.figure(figsize=(10, 5))
+plt.subplot(121)
+plot_precip_field(
+    r_radar[-1, :, :], geodata=radar_metadata, title=f"Radar observation at {date_str}"
+)
+plt.subplot(122)
+plot_precip_field(
+    r_nwp[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")
+r_radar, radar_metadata = transformer(r_radar.values, radar_metadata, threshold=0.1)
+transformer = pysteps.utils.get_method("dB")
+r_nwp, nwp_metadata = transformer(r_nwp.values, nwp_metadata, threshold=0.1)
+
+# Initial decomposition settings
+decomp_method, recomp_method = cascade.get_method("fft")
+bandpass_filter_method = "gaussian"
+M, N = r_radar.shape[1:]
+n_cascade_levels = 8
+n_models = 1  # The number of NWP models to blend with the radar rainfall nowcast
+filter_method = cascade.get_method(bandpass_filter_method)
+filter = filter_method((M, N), n_cascade_levels)
+
+# r_nwp has to be four dimentional (n_models, time, y, x).
+# If we only use one model:
+if r_nwp.ndim == 3:
+    r_nwp = r_nwp[None, :]
+
+
+################################################################################
+# Determine the velocity fields
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+oflow_method = pysteps.motion.get_method("lucaskanade")
+
+# First for the radar images
+v_radar = oflow_method(r_radar)
+
+# Then for the NWP forecast
+v_nwp = []
+# Loop through the models
+for n_model in range(r_nwp.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, r_nwp.shape[1]):
+        v_nwp_ = oflow_method(r_nwp[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)
+    v_nwp.append(_v_nwp_)
+v_nwp = np.stack(v_nwp)
+
+
+################################################################################
+# The blended forecast
+# --------------------
+
+precip_forecast = blending.steps.forecast(
+    R=r_radar,
+    R_d_models=r_nwp,
+    V=v_radar,
+    V_models=v_nwp,
+    timesteps=18,
+    timestep=timestep,
+    issuetime=date_radar,
+    n_ens_members=1,
+    n_cascade_levels=n_cascade_levels,
+    blend_nwp_members=False,
+    R_thr=radar_metadata["threshold"],
+    kmperpixel=radar_metadata["xpixelsize"] / 1000.0,
+    extrap_method="semilagrangian",
+    decomp_method="fft",
+    bandpass_filter_method="gaussian",
+    noise_method="nonparametric",
+    noise_stddev_adj="auto",
+    ar_order=2,
+    vel_pert_method=None,
+    weights_method="bps",
+    conditional=False,
+    probmatching_method="cdf",
+    mask_method="incremental",
+    callback=None,
+    return_output=True,
+    seed=None,
+    num_workers=1,
+    fft_method="numpy",
+    domain="spatial",
+    outdir_path_skill="./tmp",
+    extrap_kwargs=None,
+    filter_kwargs=None,
+    noise_kwargs=None,
+    vel_pert_kwargs=None,
+    clim_kwargs=None,
+    mask_kwargs=None,
+    measure_time=False,
+)
+
+# Transform the data back into mm/h
+precip_forecast, _ = converter(precip_forecast, radar_metadata)
+r_radar, _ = converter(r_radar, radar_metadata)
+r_nwp, _ = converter(r_nwp, 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 (t=0). Therefore, the first time
+# steps consist mostly of the extrapolation.
+# However, near the end of the forecast (t=+3h), the NWP share in the blended
+# forecast has become more important and the forecast starts to resemble the
+# NWP forecast more.
+
+# Plot the blended forecast
+plt.figure(figsize=(10, 5))
+plt.subplot(131)
+plot_precip_field(
+    precip_forecast[0, 2, :, :],
+    geodata=radar_metadata,
+    title="Blended forecast at t + 30 min",
+)
+plt.subplot(132)
+plot_precip_field(
+    precip_forecast[0, 5, :, :],
+    geodata=radar_metadata,
+    title="Blended forecast at t + 60 min",
+)
+plt.subplot(133)
+plot_precip_field(
+    precip_forecast[0, 17, :, :],
+    geodata=radar_metadata,
+    title="Blended forecast at t + 180 min",
+)
+plt.tight_layout()
+plt.show()
+
+# Plot the NWP forecast for comparison
+plt.figure(figsize=(10, 5))
+plt.subplot(131)
+plot_precip_field(
+    r_nwp[0, 3, :, :], geodata=nwp_metadata, title="NWP forecast at t + 30 min"
+)
+plt.subplot(132)
+plot_precip_field(
+    r_nwp[0, 6, :, :], geodata=nwp_metadata, title="NWP forecast at t + 60 min"
+)
+plt.subplot(133)
+plot_precip_field(
+    r_nwp[0, 18, :, :], geodata=nwp_metadata, title="NWP forecast at t + 180 min"
+)
+plt.tight_layout()
+plt.show()
+
+
+################################################################################
+# 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