Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
58 commits
Select commit Hold shift + click to select a range
c466b9c
Add importer for DWD radar products
m-rempel May 27, 2025
1e98e4f
Update DWD HDF5 Importer
m-rempel Jun 3, 2025
9d47508
Add reprojection of unstructured grids
m-rempel Jun 24, 2025
f4b1506
docs: add docstrings to function and run black
RubenImhoff Jul 1, 2025
a632102
Add suggestions of draft pull request in importer and reprojection
m-rempel Jul 18, 2025
0a78666
Merge branch 'add_dwd_radar_data' into add_enkf_combination
m-rempel Jul 21, 2025
1b063e0
Update doc strings in io/importers.py and utils/reprojection.py
m-rempel Jul 21, 2025
7d477a4
Add first version of reduced-space EnKF combination technique
m-rempel Aug 10, 2025
79427ec
Clean up EnKF code
m-rempel Aug 22, 2025
1fc2f5c
Begin documentation of EnKF code
m-rempel Aug 29, 2025
f3b9f3d
First version of reduced-space ens kalman filter code
m-rempel Sep 2, 2025
954e0f7
Move probability matching at the end of the forecast step (steps swap…
m-rempel Sep 6, 2025
7d0e2b4
fix: use isinstance instead of type
RubenImhoff Sep 9, 2025
af523b6
refactor: making a start with refactoring (to be continued)
RubenImhoff Sep 9, 2025
a6acc99
fix: fix merge issues
RubenImhoff Sep 9, 2025
3e8db51
refactor: change function docstrings and names
RubenImhoff Sep 9, 2025
7ce65f1
refactor: adjust get_precipitation_mask
RubenImhoff Sep 9, 2025
fee5871
Add suggested documentation.
m-rempel Sep 9, 2025
3d8dee0
Change names of ensemble variables in EnsembleKalmanFilter.update()
m-rempel Sep 9, 2025
5f5abfc
refactor: adjust get_lien_criterion function
RubenImhoff Sep 9, 2025
4824cc0
refactor: more refactoring of ens_kalman_filter_methods code
RubenImhoff Sep 9, 2025
fb451c0
refactor: refactor /utils/pca.py
RubenImhoff Sep 9, 2025
3dc64ea
Add use_accum_sampling_prob in combination kwargs
m-rempel Sep 9, 2025
6253aef
refactor: refactor pca_ens_kalman_filter.py
RubenImhoff Sep 9, 2025
53eed00
Merge branch 'add_reduced_space_enkf' of https://github.com/m-rempel/…
RubenImhoff Sep 9, 2025
b64208e
refactor: refactor pca_ens_kalman_filter.py - part 2
RubenImhoff Sep 9, 2025
e8acdb2
Remove a testwise print statement
m-rempel Sep 9, 2025
88ccec7
Merge branch 'master' into add_reduced_space_enkf
m-rempel Sep 9, 2025
08cd066
feat: first setup of new gallery example - to be tested in next commit
RubenImhoff Sep 9, 2025
4a130a5
Add no rain case handling
m-rempel Sep 9, 2025
32a9c6a
Add full NWP weight at the end of the combination
m-rempel Sep 9, 2025
6b84889
Refactor forecast main loop, no rain case handling as well as full NW…
m-rempel Sep 10, 2025
811af7c
refactor: minor refactor
RubenImhoff Sep 11, 2025
9a0963c
refactor: fix some issues raised by the coadacy report
RubenImhoff Sep 11, 2025
3a56108
fix: add sklearn to environment requirements
RubenImhoff Sep 11, 2025
1d160cf
fix: change examples to omit merge conflict
RubenImhoff Sep 11, 2025
fb92079
fix: set new bibliography entry to right fromat
RubenImhoff Sep 11, 2025
00e2f90
Refactor global variables, add a class for ensemble member's common m…
m-rempel Sep 11, 2025
8a9354c
Merge branch 'add_reduced_space_enkf' of github.com:m-rempel/pysteps …
m-rempel Sep 11, 2025
3441273
Merge branch 'master' into add_reduced_space_enkf
m-rempel Sep 11, 2025
4b50e5f
Run black
m-rempel Sep 11, 2025
84c898f
Resolve codacy issues
m-rempel Sep 11, 2025
0510e37
Resolve codacy issues II
m-rempel Sep 11, 2025
11ca2d3
Add importorskip('sklearn') in PCA EnKF test
m-rempel Sep 11, 2025
7b1cd8a
Add scikit-learn as an optional test dependency
m-rempel Sep 12, 2025
10ae024
minor updates
RubenImhoff Sep 12, 2025
c6d259c
Merge branch 'add_reduced_space_enkf' of https://github.com/m-rempel/…
RubenImhoff Sep 12, 2025
65fa254
Add test for PCA and resolve codacy issues
m-rempel Sep 12, 2025
89c58b9
fix: fix some minor errors
RubenImhoff Sep 12, 2025
bede95d
Merge branch 'add_reduced_space_enkf' of https://github.com/m-rempel/…
RubenImhoff Sep 12, 2025
c840e6d
Solve bug in PCA test
m-rempel Sep 12, 2025
2ff6d1b
fix: minor fixes found through testing the gallery example
RubenImhoff Sep 12, 2025
d9f2353
Merge branch 'add_reduced_space_enkf' of https://github.com/m-rempel/…
RubenImhoff Sep 12, 2025
4fa4443
feat: adjust gallery example
RubenImhoff Sep 12, 2025
91c1cbf
fix: swap indices
RubenImhoff Sep 12, 2025
fb9f841
fix: fix assertion margin in pca test
RubenImhoff Sep 12, 2025
e45ecdc
Fix bug in PCA test when n_components < n_ens_member
m-rempel Sep 12, 2025
356e924
refactor: change norain threshold in example
RubenImhoff Sep 17, 2025
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions ci/ci_test_env.yml
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ dependencies:
- PyWavelets
- pandas
- scikit-image
- scikit-learn
- rasterio
- gdal
# Test dependencies
Expand Down
2 changes: 2 additions & 0 deletions doc/source/pysteps_reference/blending.rst
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,8 @@ Implementation of blending methods for blending (ensemble) nowcasts with Numeric

.. automodule:: pysteps.blending.interface
.. automodule:: pysteps.blending.clim
.. automodule:: pysteps.blending.ens_kalman_filter_blending
.. automodule:: pysteps.blending.ens_kalman_filter_methods
.. automodule:: pysteps.blending.linear_blending
.. automodule:: pysteps.blending.skill_scores
.. automodule:: pysteps.blending.steps
Expand Down
3 changes: 2 additions & 1 deletion doc/source/pysteps_reference/utils.rst
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,8 @@ Implementation of miscellaneous utility functions.
.. automodule:: pysteps.utils.fft
.. automodule:: pysteps.utils.images
.. automodule:: pysteps.utils.interpolate
.. automodule:: pysteps.utils.pca
.. automodule:: pysteps.utils.reprojection
.. automodule:: pysteps.utils.spectral
.. automodule:: pysteps.utils.tapering
.. automodule:: pysteps.utils.transformation
.. automodule:: pysteps.utils.reprojection
11 changes: 11 additions & 0 deletions doc/source/references.bib
Original file line number Diff line number Diff line change
Expand Up @@ -388,3 +388,14 @@ @ARTICLE{Imhoff2023
YEAR = 2023,
DOI = "10.1002/qj.4461"
}

@ARTICLE{Nerini2019MWR,
title = {A {Reduced}-{Space} {Ensemble} {Kalman} {Filter} {Approach} for {Flow}-{Dependent} {Integration} of {Radar} {Extrapolation} {Nowcasts} and {NWP} {Precipitation} {Ensembles}},
volume = {147},
doi = {10.1175/MWR-D-18-0258.1},
number = {3},
journal = {Monthly Weather Review},
author = {D. Nerini and L. Foresti and D. Leuenberger and S. Robert and U. Germann},
year = {2019},
pages = {987--1006},
}
1 change: 1 addition & 0 deletions doc/source/user_guide/example_data.rst
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,7 @@ directory structure looks like this::
├── KNMI
├── OPERA
├── bom
├── dwd
├── fmi
├── mch

Expand Down
2 changes: 1 addition & 1 deletion environment.yml
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ channels:
- conda-forge
- defaults
dependencies:
- python>=3.11
- python>=3.10
- jsmin
- jsonschema
- matplotlib
Expand Down
3 changes: 2 additions & 1 deletion environment_dev.yml
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ channels:
- conda-forge
- defaults
dependencies:
- python>=3.11
- python>=3.10
- pip
- jsmin
- jsonschema
Expand All @@ -29,5 +29,6 @@ dependencies:
- pre_commit
- cartopy>=0.18
- scikit-image
- scikit-learn
- pandas
- rasterio
291 changes: 291 additions & 0 deletions examples/ens_kalman_filter_blended_forecast.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,291 @@
# -*- coding: utf-8 -*-
"""
Blended forecast
====================

This tutorial shows how to construct a blended forecast between an ensemble nowcast
and an ensemble Numerical Weather Prediction (NWP) rainfall forecast using the
Reduced-Space Ensemble Kalman filter approach in :cite:`Nerini2019MWR`. The used
datasets are from the German Weather Service (DWD).

"""

import os
from datetime import datetime, timedelta

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
import pysteps_nwp_importers


################################################################################
# Read the radar images and the NWP forecast
# ------------------------------------------
#
# First, we import a sequence of 4 images of 5-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("202506041645", "%Y%m%d%H%M")
# The last NWP forecast was issued at 16:00 - the blending tool will be able
# to find the correct lead times itself.
date_nwp = datetime.strptime("202506041600", "%Y%m%d%H%M")
radar_data_source = rcparams.data_sources["dwd"]
nwp_data_source = rcparams.data_sources["dwd_nwp"]


###############################################################################
# Load the data from the archive
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

root_path = radar_data_source["root_path"]
path_fmt = radar_data_source["path_fmt"]
fn_pattern = radar_data_source["fn_pattern"]
fn_ext = radar_data_source["fn_ext"]
importer_name = radar_data_source["importer"]
importer_kwargs = radar_data_source["importer_kwargs"]
timestep_radar = radar_data_source["timestep"]

# Find the radar files in the archive
fns = io.find_by_date(
date_radar,
root_path,
path_fmt,
fn_pattern,
fn_ext,
timestep_radar,
num_prev_files=2,
)

# Read the radar composites (which are already in mm/h)
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("dwd_nwp", "importer")
kwargs = {
"varname": "lsprate",
"grid_file_path": "./aux/grid_files/dwd/icon/R19B07/icon_grid_0047_R19B07_L.nc",
}
nwp_precip, _, nwp_metadata = nwp_importer(filename, **kwargs)


################################################################################
# Pre-processing steps
# --------------------

# Set the zerovalue and precipitation thresholds (these are fixed from DWD)
prec_thr = 0.049
zerovalue = 0.027

# Transform the zerovalue and precipitation thresholds to dBR
log_thr_prec = 10.0 * np.log10(prec_thr)
log_zerovalue = 10.0 * np.log10(zerovalue)

# Reproject the DWD ICON NWP data onto a regular grid
nwp_metadata["clon"] = nwp_precip["longitude"].values
nwp_metadata["clat"] = nwp_precip["latitude"].values
# We change the time step from the DWD NWP data to 15 min (it is actually 5 min)
# to have a longer forecast horizon available for this example, as pysteps_data
# only contains 1 hour of DWD forecast data (to minimize storage).
nwp_metadata["accutime"] = 15.0
nwp_precip = (
nwp_precip.values * 3.0
) # (to account for the change in time step from 5 to 15 min)

# Reproject ID2 data onto a regular grid
nwp_precip_rprj, nwp_metadata_rprj = (
pysteps_nwp_importers.importer_dwd_nwp.unstructured2regular(
nwp_precip, nwp_metadata, radar_metadata
)
)

# Make sure the units are in mm/h
converter = pysteps.utils.get_method("mm/h")
radar_precip, radar_metadata = converter(
radar_precip, radar_metadata
) # The radar data should already be in mm/h
nwp_precip_rprj, nwp_metadata_rprj = converter(nwp_precip_rprj, nwp_metadata_rprj)

# Threshold the data
radar_precip[radar_precip < prec_thr] = 0.0
nwp_precip_rprj[nwp_precip_rprj < prec_thr] = 0.0

# Plot the radar rainfall field and the first time step and first ensemble member
# 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}",
colorscale="STEPS-NL",
)
plt.subplot(122)
plot_precip_field(
nwp_precip_rprj[0, 0, :, :],
geodata=nwp_metadata_rprj,
title=f"NWP forecast at {date_str}",
colorscale="STEPS-NL",
)
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=prec_thr, zerovalue=log_zerovalue
)
nwp_precip_rprj, nwp_metadata_rprj = transformer(
nwp_precip_rprj, nwp_metadata_rprj, threshold=prec_thr, zerovalue=log_zerovalue
)


###############################################################################
# Determine the velocity fields
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# In contrast to the STEPS blending method, no motion field for the NWP fields
# is needed in the ensemble kalman filter blending approach.

# Estimate the motion vector field
oflow_method = pysteps.motion.get_method("lucaskanade")
velocity_radar = oflow_method(radar_precip)


################################################################################
# The blended forecast
# ~~~~~~~~~~~~~~~~~~~~

# Set the timestamps for radar_precip and nwp_precip_rprj
timestamps_radar = np.array(
sorted(
[
date_radar - timedelta(minutes=i * timestep_radar)
for i in range(len(radar_precip))
]
)
)
timestamps_nwp = np.array(
sorted(
[
date_nwp + timedelta(minutes=i * int(nwp_metadata_rprj["accutime"]))
for i in range(nwp_precip_rprj.shape[0])
]
)
)

# Set the combination kwargs
combination_kwargs = dict(
n_tapering=0, # No. of principal components of the ens. forecast for which the cov. matrix is computed
non_precip_mask=True, # Specifies whether the computation should be truncated on grid boxes where at least a minimum number of ens. members forecast precipitation.
n_ens_prec=1, # Minimum number of ens. members that forecast precip for the above-mentioned mask.
lien_criterion=True, # Specifies wheter the Lien criterion should be applied.
n_lien=10, # Minimum number of ensemble members that forecast precipitation for the Lien criterion (equals half the ens. members here)
prob_matching="iterative", # The type of probability matching used.
inflation_factor_bg=1.8, # Inflation factor of the background (NWC) covariance matrix. (this value indicates a faster convergence towards the NWP ensemble)
inflation_factor_obs=1.0, # Inflation factor of the observation (NWP) covariance matrix.
offset_bg=0.0, # Offset of the background (NWC) covariance matrix.
offset_obs=0.0, # Offset of the observation (NWP) covariance matrix.
nwp_hres_eff=14.0, # Effective horizontal resolution of the utilized NWP model (in km here).
sampling_prob_source="ensemble", # Computation method of the sampling probability for the probability matching. 'ensemble' computes this probability as the ratio between the ensemble differences.
use_accum_sampling_prob=False, # Specifies whether the current sampling probability should be used for the probability matching or a probability integrated over the previous forecast time.
)


# Call the PCA EnKF method
blending_method = blending.get_method("pca_enkf")
precip_forecast = blending_method(
obs_precip=radar_precip, # Radar data in dBR
obs_timestamps=timestamps_radar, # Radar timestamps
nwp_precip=nwp_precip_rprj, # NWP in dBR
nwp_timestamps=timestamps_nwp, # NWP timestamps
velocity=velocity_radar, # Velocity vector field
forecast_horizon=120, # Forecast length (horizon) in minutes - only a short forecast horizon due to the limited dataset length stored here.
issuetime=date_radar, # Forecast issue time as datetime object
n_ens_members=20, # No. of ensemble members
precip_mask_dilation=1, # Dilation of precipitation mask in grid boxes
n_cascade_levels=6, # No. of cascade levels
precip_thr=log_thr_prec, # Precip threshold
norain_thr=0.0005, # Minimum of 0.5% precip needed, otherwise 'zero rainfall'
num_workers=4, # No. of parallel threads
noise_stddev_adj="auto", # Standard deviation adjustment
noise_method="ssft", # SSFT as noise method
enable_combination=True, # Enable combination
noise_kwargs={"win_size": (512, 512), "win_fun": "hann", "overlap": 0.5},
extrap_kwargs={"interp_order": 3, "map_coordinates_mode": "nearest"},
combination_kwargs=combination_kwargs,
filter_kwargs={"include_mean": True},
)

# 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_rprj, nwp_metadata_rprj)


################################################################################
# Visualize the output
# ~~~~~~~~~~~~~~~~~~~~
#
# The NWP rainfall forecast has a much 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
# the more dominant contribution to the forecast and thus the forecast starts
# to resemble the NWP forecast.

fig = plt.figure(figsize=(5, 12))

leadtimes_min = [15, 30, 45, 60, 90, 120]
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_radar) - 1, :, :],
geodata=radar_metadata,
title=f"Blended +{leadtime} min",
axis="off",
colorscale="STEPS-NL",
colorbar=False,
)

# Raw NWP forecast
plt.subplot(n_leadtimes, 2, n * 2 + 2)
plot_precip_field(
nwp_precip[int(leadtime / int(nwp_metadata_rprj["accutime"])) - 1, 0, :, :],
geodata=nwp_metadata_rprj,
title=f"NWP +{leadtime} min",
axis="off",
colorscale="STEPS-NL",
colorbar=False,
)


################################################################################
# References
# ~~~~~~~~~~
#

# Nerini, D., Foresti, L., Leuenberger, D., Robert, S., Germann, U. 2019. "A
# Reduced-Space Ensemble Kalman Filter Approach for Flow-Dependent Integration
# of Radar Extrapolation Nowcasts and NWP Precipitation Ensembles." Monthly
# Weather Review 147(3): 987-1006. https://doi.org/10.1175/MWR-D-18-0258.1.
File renamed without changes.
Loading
Loading