.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/ens_kalman_filter_blended_forecast.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_ens_kalman_filter_blended_forecast.py: Ensemble-based Blending ======================= This tutorial demonstrates how to construct a blended rainfall forecast by combining an ensemble nowcast with an ensemble Numerical Weather Prediction (NWP) forecast. The method follows the Reduced-Space Ensemble Kalman Filter approach described in :cite:`Nerini2019MWR`. The procedure starts from the most recent radar observations. In the **prediction step**, a stochastic radar extrapolation technique generates short-term forecasts. In the **correction step**, these forecasts are updated using information from the latest ensemble NWP run. To make the matrix operations tractable, the Bayesian update is carried out in the subspace defined by the leading principal components—hence the term *reduced space*. The datasets used in this tutorial are provided by the German Weather Service (DWD). .. GENERATED FROM PYTHON SOURCE LINES 20-33 .. code-block:: Python 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.utils import aggregate_fields_space from pysteps.visualization import plot_precip_field import pysteps_nwp_importers .. GENERATED FROM PYTHON SOURCE LINES 34-44 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. .. GENERATED FROM PYTHON SOURCE LINES 44-54 .. code-block:: Python # 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"] .. GENERATED FROM PYTHON SOURCE LINES 55-57 Load the data from the archive ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 57-102 .. code-block:: Python 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 = nwp_data_source["importer_kwargs"] # Resolve grid_file_path relative to PYSTEPS_DATA_PATH kwargs["grid_file_path"] = os.path.join( os.environ["PYSTEPS_DATA_PATH"], kwargs["grid_file_path"] ) nwp_precip, _, nwp_metadata = nwp_importer(filename, **kwargs) # We lower the number of ens members to 10 to reduce the memory needs in the # example here. However, it is advised to have a minimum of 20 members for the # Reduced-Space Ensemble Kalman filter approach nwp_precip = nwp_precip[:, 0:10, :].astype("single") .. GENERATED FROM PYTHON SOURCE LINES 103-105 Pre-processing steps -------------------- .. GENERATED FROM PYTHON SOURCE LINES 105-186 .. code-block:: Python # 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.astype("single") * 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 ) ) nwp_precip = None # Upscale both the radar and NWP data to a twice as coarse resolution to lower # the memory needs (for this example) radar_precip, radar_metadata = aggregate_fields_space( radar_precip, radar_metadata, radar_metadata["xpixelsize"] * 4 ) nwp_precip_rprj, nwp_metadata_rprj = aggregate_fields_space( nwp_precip_rprj.astype("single"), nwp_metadata_rprj, nwp_metadata_rprj["xpixelsize"] * 4, ) # 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 ) .. image-sg:: /auto_examples/images/sphx_glr_ens_kalman_filter_blended_forecast_001.png :alt: Radar observation at 2025-06-04 16:45, NWP forecast at 2025-06-04 16:45 :srcset: /auto_examples/images/sphx_glr_ens_kalman_filter_blended_forecast_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 187-192 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. .. GENERATED FROM PYTHON SOURCE LINES 192-198 .. code-block:: Python # Estimate the motion vector field oflow_method = pysteps.motion.get_method("lucaskanade") velocity_radar = oflow_method(radar_precip) .. GENERATED FROM PYTHON SOURCE LINES 199-201 The blended forecast ~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 201-269 .. code-block:: Python # 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, # Tapering parameter: controls how many diagonals of the covariance matrix are kept (0 = no tapering) 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=5, # 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=3.0, # 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=10, # 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) .. rst-class:: sphx-glr-script-out .. code-block:: none Test for no rain cases ====================== Rain fraction is: 0.08794545454545455, while minimum fraction is 0.0005 Rain fraction is: 0.14362814814814814, while minimum fraction is 0.0005 Reduced-space ensemble Kalman filter ==================================== Inputs ------ Forecast issue time: 2025-06-04T16:45:00 Input dimensions: 300x275 Temporal resolution: 5 minutes NWP and blending inputs ----------------------- Number of (NWP) models: 10 Methods ------- Extrapolation: semilagrangian Bandpass filter: gaussian Decomposition: fft Noise generator: ssft Noise adjustment: yes EnKF implementation: masked_enkf FFT method: numpy Domain: spatial Parameters ---------- Forecast length in min: 120 Ensemble size: 10 Parallel threads: 4 Number of cascade levels: 6 Order of the AR(p) model: 1 No rain forecast: none Nowcast components initialized successfully. ************************************************ * Correlation coefficients for cascade levels: * ************************************************ ------------------------- | Level | Lag-1 | ------------------------- | 1 | 0.999399 | ------------------------- | 2 | 0.996212 | ------------------------- | 3 | 0.988787 | ------------------------- | 4 | 0.972267 | ------------------------- | 5 | 0.904330 | ------------------------- | 6 | 0.578820 | ------------------------- **************************************** * AR(p) parameters for cascade levels: * **************************************** --------------------------------------- | Level | Phi-1 | Phi-0 | --------------------------------------- | 1 | 0.999399 | 0.034674 | --------------------------------------- | 2 | 0.996212 | 0.086954 | --------------------------------------- | 3 | 0.988787 | 0.149329 | --------------------------------------- | 4 | 0.972267 | 0.233874 | --------------------------------------- | 5 | 0.904330 | 0.426835 | --------------------------------------- | 6 | 0.578820 | 0.815456 | --------------------------------------- Computing noise adjustment coefficients... noise std. dev. coeffs: [1.20439368 1.16941322 1.01525735 1.11890008 1.05759501 0.74970358] Initialize ensemble Kalman filter ================================= Non-tapered diagonals: 0 Non precip mask: True No. ens mems with precipitation: 1 Lien Criterion: True No. ens mems with precip (Lien): 5 Initialize masked ensemble Kalman filter ======================================== Iterative probability matching: True Background inflation factor: 3.0 Observation inflation factor: 1.0 Background offset: 0.0 Observation offset: 0.0 Sampling probability source: ensemble Use accum. sampling probability: False Ensure full NWP weight: True Computing combination for lead time + 0 min ...done. Computing combination for lead time + 5 min ...done. Computing combination for lead time + 10 min ...done. Computing combination for lead time + 15 min ...done. Computing combination for lead time + 20 min Sampling probability: 0.0396 ...done. Computing combination for lead time + 25 min ...done. Computing combination for lead time + 30 min ...done. Computing combination for lead time + 35 min Sampling probability: 0.0573 ...done. Computing combination for lead time + 40 min ...done. Computing combination for lead time + 45 min ...done. Computing combination for lead time + 50 min Sampling probability: 0.1833 ...done. Computing combination for lead time + 55 min ...done. Computing combination for lead time + 60 min ...done. Computing combination for lead time + 65 min Sampling probability: 0.8115 ...done. Computing combination for lead time + 70 min ...done. Computing combination for lead time + 75 min ...done. Computing combination for lead time + 80 min Sampling probability: 0.5703 ...done. Computing combination for lead time + 85 min ...done. Computing combination for lead time + 90 min ...done. Computing combination for lead time + 95 min Sampling probability: 0.8174 ...done. Computing combination for lead time + 100 min ...done. Computing combination for lead time + 105 min ...done. Computing combination for lead time + 110 min Sampling probability: 0.8645 ...done. Computing combination for lead time + 115 min ...done. Computing combination for lead time + 120 min ...done. .. GENERATED FROM PYTHON SOURCE LINES 270-279 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. .. GENERATED FROM PYTHON SOURCE LINES 279-308 .. code-block:: Python 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, ) .. image-sg:: /auto_examples/images/sphx_glr_ens_kalman_filter_blended_forecast_002.png :alt: Blended +15 min, NWP +15 min, Blended +30 min, NWP +30 min, Blended +45 min, NWP +45 min, Blended +60 min, NWP +60 min, Blended +90 min, NWP +90 min, Blended +120 min, NWP +120 min :srcset: /auto_examples/images/sphx_glr_ens_kalman_filter_blended_forecast_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 309-312 References ~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 312-317 .. code-block:: Python # 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. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 14.869 seconds) .. _sphx_glr_download_auto_examples_ens_kalman_filter_blended_forecast.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ens_kalman_filter_blended_forecast.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: ens_kalman_filter_blended_forecast.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: ens_kalman_filter_blended_forecast.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_