From b6dd04b467c88e92d3b03ca037270bcb482e582e Mon Sep 17 00:00:00 2001 From: Ruben Imhoff Date: Tue, 18 Nov 2025 08:47:35 +0100 Subject: [PATCH 01/14] feat: add smooth_radar_mask functionality from steps blending to pca_ens_kalman_filter blending --- pysteps/blending/pca_ens_kalman_filter.py | 54 +++++++++++++++++++ .../test_blending_pca_ens_kalman_filter.py | 37 +++++++------ 2 files changed, 75 insertions(+), 16 deletions(-) diff --git a/pysteps/blending/pca_ens_kalman_filter.py b/pysteps/blending/pca_ens_kalman_filter.py index e6aa076d1..058fe5189 100644 --- a/pysteps/blending/pca_ens_kalman_filter.py +++ b/pysteps/blending/pca_ens_kalman_filter.py @@ -99,6 +99,15 @@ class EnKFCombinationConfig: precip_mask_dilation: int Number of grid boxes by which the precipitation mask should be extended per timestep. + smooth_radar_mask_range: int, Default is 0. + Method to smooth the transition between the radar-NWP-noise blend and the NWP-noise + blend near the edge of the radar domain (radar mask), where the radar data is either + not present anymore or is not reliable. If set to 0 (grid cells), this generates a + normal forecast without smoothing. To create a smooth mask, this range should be a + positive value, representing a buffer band of a number of pixels by which the mask + is cropped and smoothed. The smooth radar mask removes the hard edges between NWP + and radar in the final blended product. Typically, a value between 50 and 100 km + can be used. 80 km generally gives good results. extrapolation_method: str Name of the extrapolation method to use. See the documentation of :py:mod:`pysteps.extrapolation.interface`. @@ -186,6 +195,7 @@ class EnKFCombinationConfig: precip_threshold: float | None norain_threshold: float precip_mask_dilation: int + smooth_radar_mask_range: int extrapolation_method: str decomposition_method: str bandpass_filter_method: str @@ -686,6 +696,39 @@ def run_forecast_step(self, nwp, is_correction_timestep=False): # Extrapolate the precipitation field onto the position of the current timestep. self.__advect() + # The extrapolation components are NaN outside the advected + # radar domain. This results in NaN values in the blended + # forecast outside the radar domain. Therefore, fill these + # areas with NWP, if requested. + if self.__forecast_state.config.smooth_radar_mask_range != 0: + nan_indices = np.isnan( + self.__forecast_state.nwc_prediction[self.__ens_member] + ) + # Compute the smooth dilated mask + new_mask = blending.utils.compute_smooth_dilated_mask( + nan_indices, + max_padding_size_in_px=self.__forecast_state.config.smooth_radar_mask_range, + ) + + # Ensure mask values are between 0 and 1 + mask_model = np.clip(new_mask, 0, 1) + mask_radar = np.clip(1 - new_mask, 0, 1) + + # Handle NaNs in precip_forecast_new and precip_forecast_new_mod_only by setting NaNs to 0 in the blending step + nwp_temp = np.nan_to_num(nwp, nan=0) + nwc_temp = np.nan_to_num( + self.__forecast_state.nwc_prediction[self.__ens_member], nan=0 + ) + + # Perform the blending of radar and model inside the radar domain using a weighted combination + self.__forecast_state.nwc_prediction[self.__ens_member] = np.nansum( + [ + mask_model * nwp_temp, + mask_radar * nwc_temp, + ], + axis=0, + ) + # Create the resulting precipitation field and set no data area. In future, when # transformation between linear and logarithmic scale will be necessary, it will be # implemented in this function. @@ -1471,6 +1514,7 @@ def forecast( issuetime, n_ens_members, precip_mask_dilation=1, + smooth_radar_mask_range=0, n_cascade_levels=6, precip_thr=-10.0, norain_thr=0.01, @@ -1529,6 +1573,15 @@ def forecast( precip_mask_dilation: int Range by which the precipitation mask within the forecast step should be extended per time step. Defaults to 1. + smooth_radar_mask_range: int, Default is 0. + Method to smooth the transition between the radar-NWP-noise blend and the NWP-noise + blend near the edge of the radar domain (radar mask), where the radar data is either + not present anymore or is not reliable. If set to 0 (grid cells), this generates a + normal forecast without smoothing. To create a smooth mask, this range should be a + positive value, representing a buffer band of a number of pixels by which the mask + is cropped and smoothed. The smooth radar mask removes the hard edges between NWP + and radar in the final blended product. Typically, a value between 50 and 100 km + can be used. 80 km generally gives good results. n_cascade_levels: int, optional The number of cascade levels to use. Defaults to 6, see issue #385 on GitHub. precip_thr: float, optional @@ -1646,6 +1699,7 @@ def forecast( precip_threshold=precip_thr, norain_threshold=norain_thr, precip_mask_dilation=precip_mask_dilation, + smooth_radar_mask_range=smooth_radar_mask_range, extrapolation_method=extrap_method, decomposition_method=decomp_method, bandpass_filter_method=bandpass_filter_method, diff --git a/pysteps/tests/test_blending_pca_ens_kalman_filter.py b/pysteps/tests/test_blending_pca_ens_kalman_filter.py index 65df77917..74216b0cb 100644 --- a/pysteps/tests/test_blending_pca_ens_kalman_filter.py +++ b/pysteps/tests/test_blending_pca_ens_kalman_filter.py @@ -10,37 +10,39 @@ # fmt: off pca_enkf_arg_values = [ # Standard setting - (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False), + (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), + # Smooth radar mask + (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,20), # Coarser NWP temporal resolution - (20,30,0,-60,False,False,5,15,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False), + (20,30,0,-60,False,False,5,15,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), # Coarser Obs temporal resolution - (20,30,0,-60,False,False,10,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False), + (20,30,0,-60,False,False,10,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), # Larger shift of the NWP init - (20,30,0,-30,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False), + (20,30,0,-30,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), # Zero rain case in observation - (20,30,0,-60,True,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False), + (20,30,0,-60,True,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), # Zero rain case in NWP - (20,30,0,-60,False,True,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False), + (20,30,0,-60,False,True,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), # Zero rain in both - (20,30,0,-60,True,True,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False), + (20,30,0,-60,True,True,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), # Accumulated sampling probability - (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",True,False), + (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",True,False,0), # Use full NWP weight - (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,True), + (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,True,0), # Both - (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",True,True), + (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",True,True,0), # Explained variance as sampling probability source - (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"explained_var",False,False), + (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"explained_var",False,False,0), # No combination - (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",False,None,1.0,"ensemble",False,False), + (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",False,None,1.0,"ensemble",False,False,0), # Standard deviation adjustment - (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,"auto",1.0,"ensemble",False,False), + (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,"auto",1.0,"ensemble",False,False,0), # Other number of ensemble members - (10,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False), + (10,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), # Other forecast length - (20,35,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False), + (20,35,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), # Other noise method - (20,30,0,-60,False,False,5,5,0.05,0.01,"nonparametric","masked_enkf",True,None,1.0,"ensemble",False,False),] + (20,30,0,-60,False,False,5,5,0.05,0.01,"nonparametric","masked_enkf",True,None,1.0,"ensemble",False,False,0),] # fmt: on pca_enkf_arg_names = ( @@ -62,6 +64,7 @@ "sampling_prob_source", "use_accum_sampling_prob", "ensure_full_nwp_weight", + "smooth_radar_mask_range", ) @@ -85,6 +88,7 @@ def test_pca_enkf_combination( sampling_prob_source, use_accum_sampling_prob, ensure_full_nwp_weight, + smooth_radar_mask_range, ): pytest.importorskip("sklearn") @@ -221,6 +225,7 @@ def test_pca_enkf_combination( issuetime=forecast_init, n_ens_members=n_ens_members, precip_mask_dilation=1, + smooth_radar_mask_range=smooth_radar_mask_range, n_cascade_levels=6, precip_thr=metadata["threshold"], norain_thr=norain_thr, From da9eb61aac82925e883f1fde5b91d05549abf4b9 Mon Sep 17 00:00:00 2001 From: Ruben Imhoff Date: Wed, 19 Nov 2025 18:45:22 +0100 Subject: [PATCH 02/14] fix: ensure no nans are left after forecast step and adjust smooth radar mask --- pysteps/blending/pca_ens_kalman_filter.py | 28 ++++++++++++++++------- 1 file changed, 20 insertions(+), 8 deletions(-) diff --git a/pysteps/blending/pca_ens_kalman_filter.py b/pysteps/blending/pca_ens_kalman_filter.py index 058fe5189..6e0f47843 100644 --- a/pysteps/blending/pca_ens_kalman_filter.py +++ b/pysteps/blending/pca_ens_kalman_filter.py @@ -694,19 +694,28 @@ def run_forecast_step(self, nwp, is_correction_timestep=False): self.__probability_matching() # Extrapolate the precipitation field onto the position of the current timestep. + # If smooth_radar_mask_range is not zero, ensure the extrapolation kwargs use + # a constant value instead of "nearest" for the coordinate mapping, otherwise + # there are possibly no nans in the domain. + if self.__forecast_state.config.smooth_radar_mask_range != 0: + self.__forecast_state.params.extrapolation_kwargs[ + "map_coordinates_mode" + ] = "constant" self.__advect() # The extrapolation components are NaN outside the advected # radar domain. This results in NaN values in the blended # forecast outside the radar domain. Therefore, fill these # areas with NWP, if requested. + nan_mask = np.isnan(self.__forecast_state.nwc_prediction[self.__ens_member]) + self.__forecast_state.nwc_prediction[self.__ens_member][nan_mask] = nwp + # For a smoother transition at the edge, we can slowly dilute the nowcast + # component into NWP at the edges if self.__forecast_state.config.smooth_radar_mask_range != 0: - nan_indices = np.isnan( - self.__forecast_state.nwc_prediction[self.__ens_member] - ) + nan_mask = np.isnan(self.__forecast_state.nwc_prediction[self.__ens_member]) # Compute the smooth dilated mask new_mask = blending.utils.compute_smooth_dilated_mask( - nan_indices, + nan_mask, max_padding_size_in_px=self.__forecast_state.config.smooth_radar_mask_range, ) @@ -721,7 +730,7 @@ def run_forecast_step(self, nwp, is_correction_timestep=False): ) # Perform the blending of radar and model inside the radar domain using a weighted combination - self.__forecast_state.nwc_prediction[self.__ens_member] = np.nansum( + self.__forecast_state.nwc_prediction_btf[self.__ens_member] = np.nansum( [ mask_model * nwp_temp, mask_radar * nwc_temp, @@ -732,13 +741,16 @@ def run_forecast_step(self, nwp, is_correction_timestep=False): # Create the resulting precipitation field and set no data area. In future, when # transformation between linear and logarithmic scale will be necessary, it will be # implemented in this function. + # TODO: once this transformation is needed, adjust the smoothed transition between + # radar mask and NWP as performed at the end of the run_forecast_step function. def backtransform(self): # Set the resulting field as shallow copy of the field that is used # continuously for forecast computation. - self.__forecast_state.nwc_prediction_btf[self.__ens_member] = ( - self.__forecast_state.nwc_prediction[self.__ens_member] - ) + if self.__forecast_state.config.smooth_radar_mask_range == 0: + self.__forecast_state.nwc_prediction_btf[self.__ens_member] = ( + self.__forecast_state.nwc_prediction[self.__ens_member] + ) # Set no data area self.__set_no_data() From 90888685f89c6f8d7e639caba215c0e112a8567f Mon Sep 17 00:00:00 2001 From: Ruben Imhoff <31476760+RubenImhoff@users.noreply.github.com> Date: Thu, 20 Nov 2025 11:08:02 +0100 Subject: [PATCH 03/14] fix: indices when filling with NWP Co-authored-by: mats-knmi <145579783+mats-knmi@users.noreply.github.com> --- pysteps/blending/pca_ens_kalman_filter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pysteps/blending/pca_ens_kalman_filter.py b/pysteps/blending/pca_ens_kalman_filter.py index 6e0f47843..dfdf91d48 100644 --- a/pysteps/blending/pca_ens_kalman_filter.py +++ b/pysteps/blending/pca_ens_kalman_filter.py @@ -708,7 +708,7 @@ def run_forecast_step(self, nwp, is_correction_timestep=False): # forecast outside the radar domain. Therefore, fill these # areas with NWP, if requested. nan_mask = np.isnan(self.__forecast_state.nwc_prediction[self.__ens_member]) - self.__forecast_state.nwc_prediction[self.__ens_member][nan_mask] = nwp + self.__forecast_state.nwc_prediction[self.__ens_member][nan_mask] = nwp[nan_mask] # For a smoother transition at the edge, we can slowly dilute the nowcast # component into NWP at the edges if self.__forecast_state.config.smooth_radar_mask_range != 0: From 120d8502cc77c4d8d60877eee5173b4ab8a1840a Mon Sep 17 00:00:00 2001 From: Ruben Imhoff <31476760+RubenImhoff@users.noreply.github.com> Date: Thu, 20 Nov 2025 11:09:09 +0100 Subject: [PATCH 04/14] fix: fill nans in mask Co-authored-by: mats-knmi <145579783+mats-knmi@users.noreply.github.com> --- pysteps/blending/pca_ens_kalman_filter.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pysteps/blending/pca_ens_kalman_filter.py b/pysteps/blending/pca_ens_kalman_filter.py index dfdf91d48..a38ca2b7d 100644 --- a/pysteps/blending/pca_ens_kalman_filter.py +++ b/pysteps/blending/pca_ens_kalman_filter.py @@ -718,6 +718,7 @@ def run_forecast_step(self, nwp, is_correction_timestep=False): nan_mask, max_padding_size_in_px=self.__forecast_state.config.smooth_radar_mask_range, ) + new_mask = np.nan_to_num(new_mask, nan=0) # Ensure mask values are between 0 and 1 mask_model = np.clip(new_mask, 0, 1) From 4ab35b5449e9e29bf81be2a5a81da63822695100 Mon Sep 17 00:00:00 2001 From: m-rempel Date: Thu, 20 Nov 2025 10:43:24 +0000 Subject: [PATCH 05/14] fix: correct mask for nan values outside radar domain --- pysteps/blending/pca_ens_kalman_filter.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/pysteps/blending/pca_ens_kalman_filter.py b/pysteps/blending/pca_ens_kalman_filter.py index a38ca2b7d..cdefbb586 100644 --- a/pysteps/blending/pca_ens_kalman_filter.py +++ b/pysteps/blending/pca_ens_kalman_filter.py @@ -708,14 +708,16 @@ def run_forecast_step(self, nwp, is_correction_timestep=False): # forecast outside the radar domain. Therefore, fill these # areas with NWP, if requested. nan_mask = np.isnan(self.__forecast_state.nwc_prediction[self.__ens_member]) - self.__forecast_state.nwc_prediction[self.__ens_member][nan_mask] = nwp[nan_mask] + self.__forecast_state.nwc_prediction[self.__ens_member][nan_mask] = nwp[ + nan_mask + ] # For a smoother transition at the edge, we can slowly dilute the nowcast # component into NWP at the edges if self.__forecast_state.config.smooth_radar_mask_range != 0: - nan_mask = np.isnan(self.__forecast_state.nwc_prediction[self.__ens_member]) + # nan_mask = np.isnan(self.__forecast_state.nwc_prediction[self.__ens_member]) # Compute the smooth dilated mask new_mask = blending.utils.compute_smooth_dilated_mask( - nan_mask, + self.__forecast_state.params.domain_mask, max_padding_size_in_px=self.__forecast_state.config.smooth_radar_mask_range, ) new_mask = np.nan_to_num(new_mask, nan=0) From 702e43e5aa89227ceb181c9ed50caa9837c58420 Mon Sep 17 00:00:00 2001 From: m-rempel Date: Thu, 20 Nov 2025 10:45:21 +0000 Subject: [PATCH 06/14] fix: avoid overwriting of NWP values in backtransform() --- pysteps/blending/pca_ens_kalman_filter.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pysteps/blending/pca_ens_kalman_filter.py b/pysteps/blending/pca_ens_kalman_filter.py index cdefbb586..3ec9752ef 100644 --- a/pysteps/blending/pca_ens_kalman_filter.py +++ b/pysteps/blending/pca_ens_kalman_filter.py @@ -755,8 +755,8 @@ def backtransform(self): self.__forecast_state.nwc_prediction[self.__ens_member] ) - # Set no data area - self.__set_no_data() + # Set no data area + self.__set_no_data() # Call spatial decomposition function and compute an adjusted standard deviation of # each spatial scale at timesteps where NWP information is incorporated. From 0f8b510d6b801c89279d602546ccf912e84a3bb3 Mon Sep 17 00:00:00 2001 From: m-rempel Date: Thu, 20 Nov 2025 18:58:16 +0000 Subject: [PATCH 07/14] fix: working smooth radar mask for timesteps > 2 --- pysteps/blending/pca_ens_kalman_filter.py | 78 ++++++++++++++--------- 1 file changed, 48 insertions(+), 30 deletions(-) diff --git a/pysteps/blending/pca_ens_kalman_filter.py b/pysteps/blending/pca_ens_kalman_filter.py index 3ec9752ef..20475d573 100644 --- a/pysteps/blending/pca_ens_kalman_filter.py +++ b/pysteps/blending/pca_ens_kalman_filter.py @@ -701,45 +701,26 @@ def run_forecast_step(self, nwp, is_correction_timestep=False): self.__forecast_state.params.extrapolation_kwargs[ "map_coordinates_mode" ] = "constant" + self.__nwp_for_filling = nwp self.__advect() # The extrapolation components are NaN outside the advected # radar domain. This results in NaN values in the blended # forecast outside the radar domain. Therefore, fill these - # areas with NWP, if requested. + # areas with the defined minimum value, if requested. nan_mask = np.isnan(self.__forecast_state.nwc_prediction[self.__ens_member]) - self.__forecast_state.nwc_prediction[self.__ens_member][nan_mask] = nwp[ - nan_mask - ] - # For a smoother transition at the edge, we can slowly dilute the nowcast - # component into NWP at the edges - if self.__forecast_state.config.smooth_radar_mask_range != 0: - # nan_mask = np.isnan(self.__forecast_state.nwc_prediction[self.__ens_member]) - # Compute the smooth dilated mask - new_mask = blending.utils.compute_smooth_dilated_mask( - self.__forecast_state.params.domain_mask, - max_padding_size_in_px=self.__forecast_state.config.smooth_radar_mask_range, - ) - new_mask = np.nan_to_num(new_mask, nan=0) + self.__forecast_state.nwc_prediction[self.__ens_member][nan_mask] = ( + self.__forecast_state.config.precip_threshold - 2.0 + ) - # Ensure mask values are between 0 and 1 - mask_model = np.clip(new_mask, 0, 1) - mask_radar = np.clip(1 - new_mask, 0, 1) + if self.__forecast_state.config.smooth_radar_mask_range != 0: + self.__fill_periphery_with_nwp() - # Handle NaNs in precip_forecast_new and precip_forecast_new_mod_only by setting NaNs to 0 in the blending step - nwp_temp = np.nan_to_num(nwp, nan=0) - nwc_temp = np.nan_to_num( - self.__forecast_state.nwc_prediction[self.__ens_member], nan=0 - ) + # Simple setter function to get NWP data for timestep 0. It's necessary if + # smooth_radar_mask_range is > 0 + def set_nwp_data_for_filling(self, nwp): - # Perform the blending of radar and model inside the radar domain using a weighted combination - self.__forecast_state.nwc_prediction_btf[self.__ens_member] = np.nansum( - [ - mask_model * nwp_temp, - mask_radar * nwc_temp, - ], - axis=0, - ) + self.__nwp_for_filling = nwp # Create the resulting precipitation field and set no data area. In future, when # transformation between linear and logarithmic scale will be necessary, it will be @@ -915,6 +896,38 @@ def __set_no_data(self): self.__forecast_state.params.domain_mask ] = np.nan + # Fill edge zones of the domain with NWP data if smooth_radar_mask_range is > 0 + def __fill_periphery_with_nwp(self): + + # For a smoother transition at the edge, we can slowly dilute the nowcast + # component into NWP at the edges + + # Compute the smooth dilated mask + new_mask = blending.utils.compute_smooth_dilated_mask( + self.__forecast_state.params.domain_mask, + max_padding_size_in_px=self.__forecast_state.config.smooth_radar_mask_range, + ) + new_mask = np.nan_to_num(new_mask, nan=0) + + # Ensure mask values are between 0 and 1 + mask_model = np.clip(new_mask, 0, 1) + mask_radar = np.clip(1 - new_mask, 0, 1) + + # Handle NaNs in precip_forecast_new and precip_forecast_new_mod_only by setting NaNs to 0 in the blending step + nwp_temp = np.nan_to_num(self.__nwp_for_filling, nan=0) + nwc_temp = np.nan_to_num( + self.__forecast_state.nwc_prediction[self.__ens_member], nan=0 + ) + + # Perform the blending of radar and model inside the radar domain using a weighted combination + self.__forecast_state.nwc_prediction_btf[self.__ens_member] = np.nansum( + [ + mask_model * nwp_temp, + mask_radar * nwc_temp, + ], + axis=0, + ) + class EnKFCombinationNowcaster: def __init__( @@ -1376,6 +1389,11 @@ def __integrated_nowcast_main_loop(self): self.__params.extrapolation_kwargs["return_displacement"] = True is_correction_timestep = False + [ + self.FC_Models[j].set_nwp_data_for_filling(self.__nwp_precip[j, 0]) + for j in range(self.__config.n_ens_members) + ] + for t, fc_leadtime in enumerate(self.__forecast_leadtimes): if self.__config.measure_time: starttime = time.time() From ced9c4dff89ebf78d74fba6d9826da3b383ab46b Mon Sep 17 00:00:00 2001 From: Ruben Imhoff Date: Fri, 21 Nov 2025 14:51:04 +0100 Subject: [PATCH 08/14] fix: resolve minor codacy issue --- pysteps/blending/pca_ens_kalman_filter.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/pysteps/blending/pca_ens_kalman_filter.py b/pysteps/blending/pca_ens_kalman_filter.py index 20475d573..224b854b8 100644 --- a/pysteps/blending/pca_ens_kalman_filter.py +++ b/pysteps/blending/pca_ens_kalman_filter.py @@ -1389,10 +1389,8 @@ def __integrated_nowcast_main_loop(self): self.__params.extrapolation_kwargs["return_displacement"] = True is_correction_timestep = False - [ + for j in range(self.__config.n_ens_members): self.FC_Models[j].set_nwp_data_for_filling(self.__nwp_precip[j, 0]) - for j in range(self.__config.n_ens_members) - ] for t, fc_leadtime in enumerate(self.__forecast_leadtimes): if self.__config.measure_time: From 5f2f74ea593cd51a23b15add6c460d74514f36c0 Mon Sep 17 00:00:00 2001 From: m-rempel Date: Fri, 21 Nov 2025 15:26:11 +0000 Subject: [PATCH 09/14] fix: working smooth radar mask for all timesteps --- pysteps/blending/pca_ens_kalman_filter.py | 55 +++++++++++++++-------- 1 file changed, 37 insertions(+), 18 deletions(-) diff --git a/pysteps/blending/pca_ens_kalman_filter.py b/pysteps/blending/pca_ens_kalman_filter.py index 224b854b8..9e937c713 100644 --- a/pysteps/blending/pca_ens_kalman_filter.py +++ b/pysteps/blending/pca_ens_kalman_filter.py @@ -713,15 +713,6 @@ def run_forecast_step(self, nwp, is_correction_timestep=False): self.__forecast_state.config.precip_threshold - 2.0 ) - if self.__forecast_state.config.smooth_radar_mask_range != 0: - self.__fill_periphery_with_nwp() - - # Simple setter function to get NWP data for timestep 0. It's necessary if - # smooth_radar_mask_range is > 0 - def set_nwp_data_for_filling(self, nwp): - - self.__nwp_for_filling = nwp - # Create the resulting precipitation field and set no data area. In future, when # transformation between linear and logarithmic scale will be necessary, it will be # implemented in this function. @@ -897,7 +888,7 @@ def __set_no_data(self): ] = np.nan # Fill edge zones of the domain with NWP data if smooth_radar_mask_range is > 0 - def __fill_periphery_with_nwp(self): + def fill_backtransform(self, nwp): # For a smoother transition at the edge, we can slowly dilute the nowcast # component into NWP at the edges @@ -914,7 +905,7 @@ def __fill_periphery_with_nwp(self): mask_radar = np.clip(1 - new_mask, 0, 1) # Handle NaNs in precip_forecast_new and precip_forecast_new_mod_only by setting NaNs to 0 in the blending step - nwp_temp = np.nan_to_num(self.__nwp_for_filling, nan=0) + nwp_temp = np.nan_to_num(nwp, nan=0) nwc_temp = np.nan_to_num( self.__forecast_state.nwc_prediction[self.__ens_member], nan=0 ) @@ -1263,7 +1254,7 @@ def __check_input_timestamps(self): self.__nwp_precip = np.delete( self.__nwp_precip, np.logical_or( - self.__nwp_timestamps <= self.__fc_init, + self.__nwp_timestamps < self.__fc_init, self.__nwp_timestamps > self.__fc_init + datetime.timedelta(minutes=self.__fc_period), ), @@ -1278,7 +1269,7 @@ def __check_input_timestamps(self): trunc_nwp_timestamps = ( self.__nwp_timestamps[ np.logical_and( - self.__nwp_timestamps > self.__fc_init, + self.__nwp_timestamps >= self.__fc_init, self.__nwp_timestamps <= self.__fc_init + datetime.timedelta(minutes=self.__fc_period), ) @@ -1389,9 +1380,6 @@ def __integrated_nowcast_main_loop(self): self.__params.extrapolation_kwargs["return_displacement"] = True is_correction_timestep = False - for j in range(self.__config.n_ens_members): - self.FC_Models[j].set_nwp_data_for_filling(self.__nwp_precip[j, 0]) - for t, fc_leadtime in enumerate(self.__forecast_leadtimes): if self.__config.measure_time: starttime = time.time() @@ -1438,8 +1426,39 @@ def __integrated_nowcast_main_loop(self): self.__forecast_loop(t, is_correction_timestep, is_nowcasting_timestep) # Apply back transformation - for FC_Model in self.FC_Models.values(): - FC_Model.backtransform() + if self.__config.smooth_radar_mask_range == 0: + for j, FC_Model in enumerate(self.FC_Models.values()): + FC_Model.backtransform() + else: + try: + t_fill_nwp + except NameError: + t_fill_nwp = 0 + if self.__forecast_leadtimes[t] in self.__correction_leadtimes: + t_fill_nwp = np.where( + self.__correction_leadtimes == self.__forecast_leadtimes[t] + )[0][0] + + def worker(j): + + self.FC_Models[j].fill_backtransform( + self.__nwp_precip[j, t_fill_nwp] + ) + + dask_worker_collection = [] + + if DASK_IMPORTED and self.__config.n_ens_members > 1: + for j in range(self.__config.n_ens_members): + dask_worker_collection.append(dask.delayed(worker)(j)) + dask.compute( + *dask_worker_collection, + num_workers=self.__params.num_ensemble_workers, + ) + else: + for j in range(self.__config.n_ens_members): + worker(j) + + dask_worker_collection = None self.__write_output() From 2122c4406d14e374ab8f5b027a6310e534d562ae Mon Sep 17 00:00:00 2001 From: Mats Veldhuizen Date: Wed, 26 Nov 2025 15:50:04 +0100 Subject: [PATCH 10/14] advect domain mask --- pysteps/blending/pca_ens_kalman_filter.py | 42 +++++++---------------- 1 file changed, 13 insertions(+), 29 deletions(-) diff --git a/pysteps/blending/pca_ens_kalman_filter.py b/pysteps/blending/pca_ens_kalman_filter.py index 9e937c713..b266d8b8e 100644 --- a/pysteps/blending/pca_ens_kalman_filter.py +++ b/pysteps/blending/pca_ens_kalman_filter.py @@ -49,8 +49,9 @@ forecast """ -import time + import datetime +import time from copy import deepcopy import numpy as np @@ -60,10 +61,10 @@ ) from pysteps import blending, cascade, extrapolation, noise, utils -from pysteps.nowcasts import utils as nowcast_utils -from pysteps.timeseries import autoregression, correlation from pysteps.blending.ens_kalman_filter_methods import MaskedEnKF +from pysteps.nowcasts import utils as nowcast_utils from pysteps.postprocessing import probmatching +from pysteps.timeseries import autoregression, correlation from pysteps.utils.check_norain import check_norain try: @@ -288,7 +289,6 @@ def __init__( # Initialize FFT, bandpass filters, decomposition methods, and extrapolation # method. def __initialize_nowcast_components(self): - # Initialize number of ensemble workers self.__params.num_ensemble_workers = min( self.__config.n_ens_members, @@ -431,11 +431,9 @@ def transform_to_lagrangian(precip, i): self.std_extrapolation = np.array(precip_forecast_decomp["stds"]) if self.__params.no_rain_case == "obs": - GAMMA = np.ones((self.__config.n_cascade_levels, self.__config.ar_order)) else: - # If there are values in the radar fields, compute the auto-correlations GAMMA = np.empty((self.__config.n_cascade_levels, self.__config.ar_order)) @@ -476,7 +474,6 @@ def __initialize_noise(self): self.__config.noise_method is not None and self.__params.no_rain_case != "obs" ): - # get methods for perturbations init_noise, self.__params.noise_generator = noise.get_method( self.__config.noise_method @@ -561,7 +558,6 @@ def __initialize_noise_field_pool(self): ) if self.__params.noise_generator is not None: - # Determine the noise field for each ensemble member for j in range(self.__config.n_noise_fields): epsilon = self.__params.noise_generator( @@ -596,7 +592,6 @@ def __init__( latest_obs: np.ndarray, precip_mask: np.ndarray, ): - self.config = enkf_combination_config self.params = enkf_combination_params self.noise_field_pool = noise_field_pool @@ -632,7 +627,6 @@ def __init__( sigma: np.ndarray, ens_member: int, ): - # Initialize instance variables self.__forecast_state = forecast_state self.__precip_cascades = precip_cascades @@ -670,7 +664,6 @@ def __init__( # Bundle single steps of the forecast. def run_forecast_step(self, nwp, is_correction_timestep=False): - # Decompose precipitation field. self.__decompose(is_correction_timestep) @@ -719,7 +712,6 @@ def run_forecast_step(self, nwp, is_correction_timestep=False): # TODO: once this transformation is needed, adjust the smoothed transition between # radar mask and NWP as performed at the end of the run_forecast_step function. def backtransform(self): - # Set the resulting field as shallow copy of the field that is used # continuously for forecast computation. if self.__forecast_state.config.smooth_radar_mask_range == 0: @@ -733,7 +725,6 @@ def backtransform(self): # Call spatial decomposition function and compute an adjusted standard deviation of # each spatial scale at timesteps where NWP information is incorporated. def __decompose(self, is_correction_timestep): - # Call spatial decomposition method. precip_extrap_decomp = self.__forecast_state.params.decomposition_method( self.__forecast_state.nwc_prediction[self.__ens_member], @@ -778,7 +769,6 @@ def __decompose(self, is_correction_timestep): # Call extrapolation function to extrapolate the precipitation field onto the # position of the current timestep. def __advect(self): - # Since previous displacement is the sum of displacement over all previous # timesteps, we have to compute the differences between the displacements to # get the motion vector field for one time step. @@ -796,6 +786,15 @@ def __advect(self): displacement_previous=self.__previous_displacement, **self.__forecast_state.params.extrapolation_kwargs, ) + self.__forecast_state.params.domain_mask = ( + self.__forecast_state.params.extrapolation_method( + self.__forecast_state.params.domain_mask, + self.__velocity, + [1], + interp_order=1, + outval=True, + )[0] + ) # Get the difference of the previous displacement field. self.__previous_displacement -= displacement_tmp @@ -803,7 +802,6 @@ def __advect(self): # Get a noise field out of the respective pool and iterate through the AR(1) # process. def __iterate(self): - # Get a noise field out of the noise field pool and multiply it with # precipitation mask and the standard deviation coefficients. epsilon = ( @@ -816,7 +814,6 @@ def __iterate(self): # Iterate through the AR(1) process for each cascade level. for i in range(self.__forecast_state.config.n_cascade_levels): - self.__precip_cascades[i] = autoregression.iterate_ar_model( self.__precip_cascades[i], self.__forecast_state.params.PHI[i], @@ -826,7 +823,6 @@ def __iterate(self): # Update the precipitation mask for the forecast step by incorporating areas # where the NWP model forecast precipitation. def __update_precip_mask(self, nwp): - # Get the area where the NWP ensemble member forecast precipitation above # precipitation threshold and dilate it by a configurable range. precip_mask = ( @@ -871,7 +867,6 @@ def __update_precip_mask(self, nwp): # Apply probability matching def __probability_matching(self): - # Apply probability matching self.__forecast_state.nwc_prediction[self.__ens_member] = ( probmatching.nonparam_match_empirical_cdf( @@ -882,14 +877,12 @@ def __probability_matching(self): # Set no data area in the resulting precipitation field. def __set_no_data(self): - self.__forecast_state.nwc_prediction_btf[self.__ens_member][ self.__forecast_state.params.domain_mask ] = np.nan # Fill edge zones of the domain with NWP data if smooth_radar_mask_range is > 0 def fill_backtransform(self, nwp): - # For a smoother transition at the edge, we can slowly dilute the nowcast # component into NWP at the edges @@ -1281,7 +1274,6 @@ def __check_input_timestamps(self): ) def __check_no_rain_case(self): - print("Test for no rain cases") print("======================") print("") @@ -1373,7 +1365,6 @@ def __print_forecast_info(self): print(f"No rain forecast: {self.__params.no_rain_case}") def __integrated_nowcast_main_loop(self): - if self.__config.measure_time: starttime_mainloop = time.time() @@ -1406,7 +1397,6 @@ def __integrated_nowcast_main_loop(self): # If full NWP weight is reached, set pure NWP ensemble forecast in combined # forecast output if is_full_nwp_weight: - # Set t_corr to the first available NWP data timestep and that is 0 try: t_corr @@ -1440,7 +1430,6 @@ def __integrated_nowcast_main_loop(self): )[0][0] def worker(j): - self.FC_Models[j].fill_backtransform( self.__nwp_precip[j, t_fill_nwp] ) @@ -1476,7 +1465,6 @@ def worker(j): return def __forecast_loop(self, t, is_correction_timestep, is_nowcasting_timestep): - # If the temporal resolution of the NWP data is equal to those of the # observation, the correction step can be applied after the forecast # step for the current forecast leadtime. @@ -1498,7 +1486,6 @@ def __forecast_loop(self, t, is_correction_timestep, is_nowcasting_timestep): # Run nowcasting time step if is_nowcasting_timestep: - # Set t_corr to the first available NWP data timestep and that is 0 try: t_corr @@ -1506,7 +1493,6 @@ def __forecast_loop(self, t, is_correction_timestep, is_nowcasting_timestep): t_corr = 0 def worker(j): - self.FC_Models[j].run_forecast_step( nwp=self.__nwp_precip[j, t_corr], is_correction_timestep=is_correction_timestep, @@ -1528,7 +1514,6 @@ def worker(j): dask_worker_collection = None def __write_output(self): - if ( self.__config.callback is not None and self.FS.nwc_prediction_btf.shape[1] > 0 @@ -1536,7 +1521,6 @@ def __write_output(self): self.__config.callback(self.FS.nwc_prediction_btf) if self.__config.return_output: - self.FS.final_combined_forecast.append(self.FS.nwc_prediction_btf.copy()) def __measure_time(self, label, start_time): From 07725cad3d9f61b937b61b2b41868eb5187ce110 Mon Sep 17 00:00:00 2001 From: Ruben Imhoff Date: Wed, 26 Nov 2025 15:57:05 +0100 Subject: [PATCH 11/14] fix: remove unnecessary line and only advect domain mask if smooth_radar_mask_range is larger than zero --- pysteps/blending/pca_ens_kalman_filter.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/pysteps/blending/pca_ens_kalman_filter.py b/pysteps/blending/pca_ens_kalman_filter.py index b266d8b8e..b68fd026a 100644 --- a/pysteps/blending/pca_ens_kalman_filter.py +++ b/pysteps/blending/pca_ens_kalman_filter.py @@ -694,7 +694,6 @@ def run_forecast_step(self, nwp, is_correction_timestep=False): self.__forecast_state.params.extrapolation_kwargs[ "map_coordinates_mode" ] = "constant" - self.__nwp_for_filling = nwp self.__advect() # The extrapolation components are NaN outside the advected @@ -786,15 +785,16 @@ def __advect(self): displacement_previous=self.__previous_displacement, **self.__forecast_state.params.extrapolation_kwargs, ) - self.__forecast_state.params.domain_mask = ( - self.__forecast_state.params.extrapolation_method( - self.__forecast_state.params.domain_mask, - self.__velocity, - [1], - interp_order=1, - outval=True, - )[0] - ) + if self.__forecast_state.config.smooth_radar_mask_range > 0: + self.__forecast_state.params.domain_mask = ( + self.__forecast_state.params.extrapolation_method( + self.__forecast_state.params.domain_mask, + self.__velocity, + [1], + interp_order=1, + outval=True, + )[0] + ) # Get the difference of the previous displacement field. self.__previous_displacement -= displacement_tmp From d0f5b8ab4c8ab09687f8d52215c44f7d4c10bd67 Mon Sep 17 00:00:00 2001 From: m-rempel Date: Thu, 11 Dec 2025 12:59:20 +0000 Subject: [PATCH 12/14] fix: use correct NWP timesteps when complete transition to NWP data and advect domain mask only once --- pysteps/blending/pca_ens_kalman_filter.py | 26 +++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/pysteps/blending/pca_ens_kalman_filter.py b/pysteps/blending/pca_ens_kalman_filter.py index b68fd026a..da57ee2ae 100644 --- a/pysteps/blending/pca_ens_kalman_filter.py +++ b/pysteps/blending/pca_ens_kalman_filter.py @@ -176,6 +176,8 @@ class EnKFCombinationConfig: :py:mod:`pysteps.blending.ens_kalman_filter_methods`. measure_time: bool If set to True, measure, print and return the computation time. + verbose_output: bool + If set to True, return additionally the background ensemble of the EnKF for further statistics. callback: function, optional Optional function that is called after computation of each time step of the nowcast. The function takes one argument: a three-dimensional array @@ -214,6 +216,7 @@ class EnKFCombinationConfig: noise_kwargs: dict[str, Any] = field(default_factory=dict) combination_kwargs: dict[str, Any] = field(default_factory=dict) measure_time: bool = False + verbose_output: bool = False callback: Any | None = None return_output: bool = True n_noise_fields: int = 30 @@ -609,6 +612,7 @@ def __init__( self.nwc_prediction_btf = self.nwc_prediction.copy() self.final_combined_forecast = [] + self.background_ensemble = {} return @@ -785,7 +789,10 @@ def __advect(self): displacement_previous=self.__previous_displacement, **self.__forecast_state.params.extrapolation_kwargs, ) - if self.__forecast_state.config.smooth_radar_mask_range > 0: + if ( + self.__forecast_state.config.smooth_radar_mask_range > 0 + and self.__ens_member == 0 + ): self.__forecast_state.params.domain_mask = ( self.__forecast_state.params.extrapolation_method( self.__forecast_state.params.domain_mask, @@ -1082,7 +1089,9 @@ def compute_forecast(self): self.__fc_init, self.__mainloop_time, ) - return self.FS.final_combined_forecast + if self.__config.verbose_output: + return self.FS.final_combined_forecast, self.FS.background_ensemble + return self.FC.final_combined_forecast # Else, return None return None @@ -1406,7 +1415,7 @@ def __integrated_nowcast_main_loop(self): print(f"Full NWP weight is reached for lead time + {fc_leadtime} min") if is_correction_timestep: t_corr = np.where( - self.__correction_leadtimes == self.__forecast_leadtimes[t - 1] + self.__correction_leadtimes == self.__forecast_leadtimes[t] )[0][0] self.FS.nwc_prediction = self.__nwp_precip[:, t_corr] @@ -1476,6 +1485,11 @@ def __forecast_loop(self, t, is_correction_timestep, is_nowcasting_timestep): self.__correction_leadtimes == self.__forecast_leadtimes[t - 1] )[0][0] + if self.__config.verbose_output: + self.FS.background_ensemble[self.__correction_leadtimes[t_corr]] = ( + self.FS.nwc_prediction.copy() + ) + self.FS.nwc_prediction, self.FS.fc_resampled = ( self.KalmanFilterModel.correct_step( self.FS.nwc_prediction, @@ -1571,6 +1585,7 @@ def forecast( noise_kwargs=None, combination_kwargs=None, measure_time=False, + verbose_output=False, ): """ Generate a combined nowcast ensemble by using the reduced-space ensemble Kalman @@ -1700,7 +1715,9 @@ def forecast( ensemble Kalman filter method. See the documentation of :py:mod:`pysteps.blending.ens_kalman_filter_methods`. measure_time: bool - If set to True, measure, print and return the computation time. + If set to True, measure, print and return the computation time. + verbose_output: bool + If set to True, return additionally the background ensemble of the EnKF for further statistics. Returns ------- @@ -1751,6 +1768,7 @@ def forecast( noise_kwargs=noise_kwargs, combination_kwargs=combination_kwargs, measure_time=measure_time, + verbose_output=verbose_output, callback=callback, return_output=return_output, n_noise_fields=30, From a5061d9039f01b783813c19136611d5d87d885c1 Mon Sep 17 00:00:00 2001 From: m-rempel Date: Fri, 12 Dec 2025 09:56:31 +0000 Subject: [PATCH 13/14] fix: typo --- pysteps/blending/pca_ens_kalman_filter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pysteps/blending/pca_ens_kalman_filter.py b/pysteps/blending/pca_ens_kalman_filter.py index da57ee2ae..dcab310ca 100644 --- a/pysteps/blending/pca_ens_kalman_filter.py +++ b/pysteps/blending/pca_ens_kalman_filter.py @@ -1091,7 +1091,7 @@ def compute_forecast(self): ) if self.__config.verbose_output: return self.FS.final_combined_forecast, self.FS.background_ensemble - return self.FC.final_combined_forecast + return self.FS.final_combined_forecast # Else, return None return None From 1024bdc9c9dcf66c3ea115e24e1bc47434544ab0 Mon Sep 17 00:00:00 2001 From: m-rempel Date: Fri, 12 Dec 2025 10:32:56 +0000 Subject: [PATCH 14/14] fix: extend unit test for additional output --- .../test_blending_pca_ens_kalman_filter.py | 43 +++++++++++-------- 1 file changed, 26 insertions(+), 17 deletions(-) diff --git a/pysteps/tests/test_blending_pca_ens_kalman_filter.py b/pysteps/tests/test_blending_pca_ens_kalman_filter.py index 74216b0cb..3ed29bb1a 100644 --- a/pysteps/tests/test_blending_pca_ens_kalman_filter.py +++ b/pysteps/tests/test_blending_pca_ens_kalman_filter.py @@ -10,39 +10,41 @@ # fmt: off pca_enkf_arg_values = [ # Standard setting - (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), + (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0,False), # Smooth radar mask - (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,20), + (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,20,False), # Coarser NWP temporal resolution - (20,30,0,-60,False,False,5,15,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), + (20,30,0,-60,False,False,5,15,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0,False), # Coarser Obs temporal resolution - (20,30,0,-60,False,False,10,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), + (20,30,0,-60,False,False,10,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0,False), # Larger shift of the NWP init - (20,30,0,-30,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), + (20,30,0,-30,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0,False), # Zero rain case in observation - (20,30,0,-60,True,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), + (20,30,0,-60,True,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0,False), # Zero rain case in NWP - (20,30,0,-60,False,True,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), + (20,30,0,-60,False,True,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0,False), # Zero rain in both - (20,30,0,-60,True,True,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), + (20,30,0,-60,True,True,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0,False), # Accumulated sampling probability - (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",True,False,0), + (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",True,False,0,False), # Use full NWP weight - (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,True,0), + (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,True,0,False), # Both - (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",True,True,0), + (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",True,True,0,False), # Explained variance as sampling probability source - (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"explained_var",False,False,0), + (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"explained_var",False,False,0,False), # No combination - (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",False,None,1.0,"ensemble",False,False,0), + (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",False,None,1.0,"ensemble",False,False,0,False), # Standard deviation adjustment - (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,"auto",1.0,"ensemble",False,False,0), + (20,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,"auto",1.0,"ensemble",False,False,0,False), # Other number of ensemble members - (10,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), + (10,30,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0,False), # Other forecast length - (20,35,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0), + (20,35,0,-60,False,False,5,5,0.05,0.01,"ssft","masked_enkf",True,None,1.0,"ensemble",False,False,0,False), # Other noise method - (20,30,0,-60,False,False,5,5,0.05,0.01,"nonparametric","masked_enkf",True,None,1.0,"ensemble",False,False,0),] + (20,30,0,-60,False,False,5,5,0.05,0.01,"nonparametric","masked_enkf",True,None,1.0,"ensemble",False,False,0,False), + # Verbose output + (20,30,0,-60,False,False,5,5,0.05,0.01,"nonparametric","masked_enkf",True,None,1.0,"ensemble",False,False,0,True),] # fmt: on pca_enkf_arg_names = ( @@ -65,6 +67,7 @@ "use_accum_sampling_prob", "ensure_full_nwp_weight", "smooth_radar_mask_range", + "verbose_output", ) @@ -89,6 +92,7 @@ def test_pca_enkf_combination( use_accum_sampling_prob, ensure_full_nwp_weight, smooth_radar_mask_range, + verbose_output, ): pytest.importorskip("sklearn") @@ -248,8 +252,13 @@ def test_pca_enkf_combination( noise_kwargs=None, combination_kwargs=combination_kwargs, measure_time=False, + verbose_output=verbose_output, ) + if verbose_output: + assert len(combined_forecast) == 2, "Wrong amount of output data" + combined_forecast = combined_forecast[0] + assert combined_forecast.ndim == 4, "Wrong amount of dimensions in forecast output" assert ( combined_forecast.shape[0] == n_ens_members