From bff866763eb66e21d83b7d63567e90e7cba0e4e2 Mon Sep 17 00:00:00 2001 From: Joep1999 Date: Thu, 11 Sep 2025 16:59:06 +0200 Subject: [PATCH 01/18] started to change all functions to work with 4d data --- pysteps/blending/steps.py | 103 ++++++++++++++++++++++++++++---------- 1 file changed, 77 insertions(+), 26 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 0e8f28785..dadf6f91c 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -795,7 +795,11 @@ def __check_inputs(self): "precip_models must be either a two-dimensional array containing dictionaries with decomposed model fields" "or a four-dimensional array containing the original (NWP) model forecasts" ) - + precip_nowcast_dim = self.__precip_nowcast.ndim + if precip_nowcast_dim != 4: + raise ValueError( + "precip_nowcast must be a four-dimensional array containing the externally calculated nowcast" + ) if self.__config.extrapolation_kwargs is None: self.__state.extrapolation_kwargs = dict() else: @@ -893,7 +897,7 @@ def __print_forecast_info(self): ) if self.__precip_nowcast is not None: print( - f"input dimensions pre-computed nowcast: {self.__precip_nowcast.shape[1]}x{self.__precip_nowcast.shape[2]}" + f"input dimensions pre-computed nowcast: {self.__precip_nowcast.shape[2]}x{self.__precip_nowcast.shape[3]}" ) if self.__config.kmperpixel is not None: print(f"km/pixel: {self.__config.kmperpixel}") @@ -1083,11 +1087,12 @@ def transform_to_lagrangian(precip, i): ) if self.__precip_nowcast is not None: self.__precip_nowcast = self.__precip_nowcast.copy() - # TODO: check this function for ensemble nowcasts - for i in range(self.__precip_nowcast.shape[0]): - self.__precip_nowcast[i, ~np.isfinite(self.__precip_nowcast[i, :])] = ( - np.nanmin(self.__precip_nowcast[i, :]) - ) + # TODO: check this function for ensemble nowcasts -> now made it always work with 4 dimentions + for m in range(self.__precip_nowcast.shape[0]): + for t in range(self.__precip_nowcast.shape[1]): + self.__precip_nowcast[ + m, t, ~np.isfinite(self.__precip_nowcast[m, t, :, :]) + ] = np.nanmin(self.__precip_nowcast[m, t, :, :]) # Perform the cascade decomposition for the input precip fields and, # if necessary, for the (NWP) model fields @@ -1109,25 +1114,65 @@ def transform_to_lagrangian(precip, i): # Decompose precomputed nowcasts and rearange them again into the required components # TODO: check if this also works if an ensemble nowcast is provided! if self.__precip_nowcast is not None: - precip_nowcast_decomp = [] - for i in range(self.__precip_nowcast.shape[0]): - cascades = self.__params.decomposition_method( - self.__precip_nowcast[i], - self.__params.bandpass_filter, - n_levels=self.__config.n_cascade_levels, - mask=self.__params.mask_threshold, - method="fft", - fft_method=self.__params.fft, - output_domain=self.__config.domain, - compute_stats=True, - normalize=True, - compact_output=True, - ) - - precip_nowcast_decomp.append(cascades) + # precip_nowcast_decomp = [] + # for i in range(self.__precip_nowcast.shape[0]): + # cascades = self.__params.decomposition_method( + # self.__precip_nowcast[i], + # self.__params.bandpass_filter, + # n_levels=self.__config.n_cascade_levels, + # mask=self.__params.mask_threshold, + # method="fft", + # fft_method=self.__params.fft, + # output_domain=self.__config.domain, + # compute_stats=True, + # normalize=True, + # compact_output=True, + # ) + + # precip_nowcast_decomp.append(cascades) + + # self.__state.precip_nowcast_cascades = nowcast_utils.stack_cascades( + # precip_nowcast_decomp, self.__config.n_cascade_levels + # ) + if self.__precip_nowcast.shape[0] == 1: + decomp_precip_nowcast = [ + self.__params.decomposition_method( + self.__precip_nowcast[0, t, :, :], + bp_filter=self.__params.bandpass_filter, + n_levels=self.__config.n_cascade_levels, + mask=self.__params.mask_threshold, + method="fft", + fft_method=self.__params.fft, + output_domain=self.__config.domain, + compute_stats=True, + normalize=True, + compact_output=True, + ) + ] + else: + with ThreadPool(self.__config.num_workers) as pool: + decomp_precip_nowcast = pool.map( + partial( + self.__params.decomposition_method, + bp_filter=self.__params.bandpass_filter, + fft_method=self.__params.fft, + output_domain=self.__config.domain, + normalize=True, + compute_stats=True, + compact_output=True, + ), + list(self.__precip_nowcast[:, t, :, :]), + ) - self.__state.precip_nowcast_cascades = nowcast_utils.stack_cascades( - precip_nowcast_decomp, self.__config.n_cascade_levels + self.__state.precip_nowcast_cascades = np.array( + [decomp["cascade_levels"] for decomp in decomp_precip_nowcast] + ) + # TODO: check if this is also necessary for ensemble nowcasts somwhere, otherwise remove it + self.__state.mean_nowcast_timestep = np.array( + [decomp["means"] for decomp in decomp_precip_nowcast] + ) + self.__state.std_nowcast_timestep = np.array( + [decomp["stds"] for decomp in decomp_precip_nowcast] ) # Rearrange the cascaded into a four-dimensional array of shape @@ -1729,6 +1774,12 @@ def __find_nowcast_NWP_combination(self, t): "The number of NWP model members is larger than the given number of ensemble members. n_model_members <= n_ens_members." ) + n_ens_members_provided = self.__precip_nowcast.shape[0] + if n_ens_members_provided > self.__config.n_ens_members: + raise ValueError( + "The number of nowcast ensemble members provided is larger than the given number of ensemble members requested. n_ens_members_provided <= n_ens_members." + ) + # Check if NWP models/members should be used individually, or if all of # them are blended together per nowcast ensemble member. if self.__config.blend_nwp_members: @@ -2021,7 +2072,7 @@ def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): for i in range(self.__config.n_cascade_levels): # Use a deterministic Externally computed nowcasting model worker_state.precip_cascades[j][i] = ( - self.__state.precip_nowcast_cascades[i][t] + self.__state.precip_nowcast_cascades[j][i][t] ) # Follow the 'standard' STEPS blending approach as described in :cite:`Imhoff2023` From 62ab9c114547a85d51cd3b8d700ae032c7f51958 Mon Sep 17 00:00:00 2001 From: Joep1999 Date: Fri, 12 Sep 2025 17:38:24 +0200 Subject: [PATCH 02/18] inbetween commit, things not working currently --- pysteps/blending/steps.py | 183 +++++++++++++++++++++++++++++++++----- 1 file changed, 160 insertions(+), 23 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index dadf6f91c..00fe1cca0 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -1135,7 +1135,7 @@ def transform_to_lagrangian(precip, i): # precip_nowcast_decomp, self.__config.n_cascade_levels # ) if self.__precip_nowcast.shape[0] == 1: - decomp_precip_nowcast = [ + precip_nowcast_decomp_timestep = [ self.__params.decomposition_method( self.__precip_nowcast[0, t, :, :], bp_filter=self.__params.bandpass_filter, @@ -1151,7 +1151,7 @@ def transform_to_lagrangian(precip, i): ] else: with ThreadPool(self.__config.num_workers) as pool: - decomp_precip_nowcast = pool.map( + precip_nowcast_decomp_timestep = pool.map( partial( self.__params.decomposition_method, bp_filter=self.__params.bandpass_filter, @@ -1164,16 +1164,20 @@ def transform_to_lagrangian(precip, i): list(self.__precip_nowcast[:, t, :, :]), ) - self.__state.precip_nowcast_cascades = np.array( - [decomp["cascade_levels"] for decomp in decomp_precip_nowcast] + precip_nowcast_cascades = np.array( + [decomp["cascade_levels"] for decomp in precip_nowcast_decomp_timestep] ) - # TODO: check if this is also necessary for ensemble nowcasts somwhere, otherwise remove it - self.__state.mean_nowcast_timestep = np.array( - [decomp["means"] for decomp in decomp_precip_nowcast] + + self.__state.precip_nowcast_decomp = nowcast_utils.stack_cascades( + precip_nowcast_cascades, self.__config.n_cascade_levels ) - self.__state.std_nowcast_timestep = np.array( - [decomp["stds"] for decomp in decomp_precip_nowcast] + + # TODO: check if this is also necessary for ensemble nowcasts somwhere, otherwise remove it + precip_nowcast_cascades = precip_nowcast_cascades[-1] + self.__state.mean_nowcast = np.array( + np.array(precip_nowcast_cascades["means"]) ) + self.__state.std_nowcast = np.array(precip_nowcast_cascades["stds"]) # Rearrange the cascaded into a four-dimensional array of shape # (n_cascade_levels,ar_order+1,m,n) for the autoregressive model @@ -1204,13 +1208,15 @@ def transform_to_lagrangian(precip, i): self.__precip_models = np.stack(temp_precip_models) - # Check for zero input fields in the radar and NWP data. + # Check for zero input fields in the radar, nowcast and NWP data. self.__params.zero_precip_radar = check_norain( self.__precip, self.__params.precip_threshold, self.__config.norain_threshold, self.__params.noise_kwargs["win_fun"], ) + + # TODO This variable is currently not used anywhere. Should we calculate it? # The norain fraction threshold used for nwp is the default value of 0.0, # since nwp does not suffer from clutter. self.__params.zero_precip_model_fields = check_norain( @@ -1764,6 +1770,10 @@ def __find_nowcast_NWP_combination(self, t): With the way it is implemented at this moment: n_ens_members of the output equals the maximum number of (ensemble) members in the input (either the nowcasts or NWP). """ + self.__state.precip_nowcast_timestep = self.__precip_nowcast[:, t, :, :].astype( + np.float64, copy=False + ) + self.__state.velocity_models_timestep = self.__velocity_models[ :, t, :, :, : ].astype(np.float64, copy=False) @@ -1785,6 +1795,132 @@ def __find_nowcast_NWP_combination(self, t): if self.__config.blend_nwp_members: self.__state.mapping_list_NWP_member_to_ensemble_member = None + elif self.__config.nowcasting_method == "external_nowcast": + n_ens_members_max = max(n_ens_members_provided, n_model_members) + n_ens_members_min = min(n_ens_members_provided, n_model_members) + print("Number of nowcast ensembles:", n_ens_members_provided) + print("Number of nwp ensembles:", n_model_members) + # Also make a list of the model index numbers. These indices are needed + # for indexing the right climatological skill file when pysteps calculates + # the blended forecast in parallel. + if n_model_members > 1: + self.__state.mapping_list_NWP_member_to_ensemble_member = np.arange( + n_model_members + ) + else: + self.__state.mapping_list_NWP_member_to_ensemble_member = [0] + + # Now, repeat the nowcast ensemble members or the nwp models/members until + # it has the same amount of members as n_ens_members_max. For instance, if + # you have 10 ensemble nowcasts members and 3 NWP members, the output will + # be an ensemble of 10 members. Hence, the three NWP members are blended + # with the first three members of the nowcast (member one with member one, + # two with two, etc.), subsequently, the same NWP members are blended with + # the next three members (NWP member one with member 4, NWP member 2 with + # member 5, etc.), until 10 is reached. + if n_ens_members_min != n_ens_members_max: + if n_model_members == 1: + print("Repeating the NWP model for all ensemble members") + self.__state.precip_models_cascades_timestep = np.repeat( + self.__state.precip_models_cascades_timestep, + n_ens_members_max, + axis=0, + ) + self.__state.mean_models_timestep = np.repeat( + self.__state.mean_models_timestep, n_ens_members_max, axis=0 + ) + self.__state.std_models_timestep = np.repeat( + self.__state.std_models_timestep, n_ens_members_max, axis=0 + ) + self.__state.velocity_models_timestep = np.repeat( + self.__state.velocity_models_timestep, n_ens_members_max, axis=0 + ) + # For the prob. matching + self.__state.precip_models_timestep = np.repeat( + self.__state.precip_models_timestep, n_ens_members_max, axis=0 + ) + # Finally, for the model indices + self.__state.mapping_list_NWP_member_to_ensemble_member = np.repeat( + self.__state.mapping_list_NWP_member_to_ensemble_member, + n_ens_members_max, + axis=0, + ) + elif n_ens_members_provided == 1: + print("Repeating the nowcast for all ensemble members") + self.__state.precip_nowcast_cascades = np.repeat( + self.__state.precip_nowcast_cascades, + n_ens_members_max, + axis=0, + ) + self.__state.mean_nowcast = np.repeat( + self.__state.mean_nowcast, n_ens_members_max, axis=0 + ) + self.__state.std_nowcast_timestep = np.repeat( + self.__state.std_nowcast, n_ens_members_max, axis=0 + ) + self.__state.velocity_nowcast = np.repeat( + self.__state.velocity_nowcast, n_ens_members_max, axis=0 + ) + # For the prob. matching + self.__state.precip_nowcast_timestep = np.repeat( + self.__state.precip_nowcast_timestep, n_ens_members_max, axis=0 + ) + + elif n_model_members == n_ens_members_min: + print("Repeating the NWP model for all ensemble members") + repeats = [ + (n_ens_members_max + i) // n_ens_members_min + for i in range(n_ens_members_min) + ] + self.__state.precip_models_cascades_timestep = np.repeat( + self.__state.precip_models_cascades_timestep, + repeats, + axis=0, + ) + self.__state.mean_models_timestep = np.repeat( + self.__state.mean_models_timestep, repeats, axis=0 + ) + self.__state.std_models_timestep = np.repeat( + self.__state.std_models_timestep, repeats, axis=0 + ) + self.__state.velocity_models_timestep = np.repeat( + self.__state.velocity_models_timestep, repeats, axis=0 + ) + # For the prob. matching + self.__state.precip_models_timestep = np.repeat( + self.__state.precip_models_timestep, repeats, axis=0 + ) + # Finally, for the model indices + self.__state.mapping_list_NWP_member_to_ensemble_member = np.repeat( + self.__state.mapping_list_NWP_member_to_ensemble_member, + repeats, + axis=0, + ) + elif n_ens_members_provided == n_ens_members_min: + print("Repeating the nowcast for all ensemble members") + repeats = [ + (n_ens_members_max + i) // n_ens_members_min + for i in range(n_ens_members_min) + ] + self.__state.precip_nowcast_cascades = np.repeat( + self.__state.precip_nowcast_cascades, + repeats, + axis=0, + ) + self.__state.mean_nowcast = np.repeat( + self.__state.mean_nowcast, repeats, axis=0 + ) + self.__state.std_nowcast = np.repeat( + self.__state.std_nowcast, repeats, axis=0 + ) + self.__state.velocity_nowcast = np.repeat( + self.__state.velocity_nowcast, repeats, axis=0 + ) + # For the prob. matching + self.__state.precip_nowcast_timestep = np.repeat( + self.__state.precip_nowcast_timestep, repeats, axis=0 + ) + else: # Start with determining the maximum and mimimum number of members/models # in both input products @@ -2065,15 +2201,15 @@ def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): # If nowcast method seleced is external_nowcast, n_ens members has to be 1 if self.__config.nowcasting_method == "external_nowcast": print("Using nowcasting method:", self.__config.nowcasting_method) - if self.__config.n_ens_members != 1: - raise ValueError( - "Currently only 1 ensemble member supported when using (deterministic) external_nowcast as input." - ) + # if self.__config.n_ens_members != 1: + # raise ValueError( + # "Currently only 1 ensemble member supported when using (deterministic) external_nowcast as input." + # ) for i in range(self.__config.n_cascade_levels): # Use a deterministic Externally computed nowcasting model - worker_state.precip_cascades[j][i] = ( - self.__state.precip_nowcast_cascades[j][i][t] - ) + worker_state.precip_cascades[j][i] = self.__state.precip_nowcast_decomp[ + j + ][i][t] # Follow the 'standard' STEPS blending approach as described in :cite:`Imhoff2023` elif self.__config.nowcasting_method == "steps": @@ -2142,8 +2278,9 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( worker_state.precip_extrapolated_probability_matching = [] if ( - self.__config.nowcasting_method == "external_nowcast" - and self.__config.n_ens_members == 1 + self.__config.nowcasting_method + == "external_nowcast" + # and self.__config.n_ens_members == 1 ): for i in range(self.__config.n_cascade_levels): precip_extrapolated_decomp = worker_state.precip_cascades[j][i][-1, :] @@ -2541,14 +2678,14 @@ def __blend_cascades(self, t_sub, j, worker_state): ) # [(extr_field, n_model_fields), n_cascade_levels, ...] means_stacked = np.concatenate( ( - worker_state.mean_extrapolation[None, :], + worker_state.mean_nowcast[None, :], worker_state.mean_models_timestep, ), axis=0, ) sigmas_stacked = np.concatenate( ( - worker_state.std_extrapolation[None, :], + worker_state.std_nowcast[None, :], worker_state.std_models_timestep, ), axis=0, @@ -2597,14 +2734,14 @@ def __blend_cascades(self, t_sub, j, worker_state): ) # [(extr_field, n_model_fields), n_cascade_levels, ...] means_stacked = np.concatenate( ( - worker_state.mean_extrapolation[None, :], + worker_state.mean_nowcast[None, :], worker_state.mean_models_timestep, ), axis=0, ) sigmas_stacked = np.concatenate( ( - worker_state.std_extrapolation[None, :], + worker_state.std_nowcast[None, :], worker_state.std_models_timestep, ), axis=0, From 099db2c2bc4ca3a38e3dc6e1a4ac2319f25bd407 Mon Sep 17 00:00:00 2001 From: Joep1999 Date: Mon, 15 Sep 2025 17:36:28 +0200 Subject: [PATCH 03/18] ensembles working, but blending does not look good --- pysteps/blending/steps.py | 132 +++++++++++++++++++++++--------------- 1 file changed, 80 insertions(+), 52 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 00fe1cca0..95faeb8dc 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -352,7 +352,6 @@ class StepsBlendingParams: class StepsBlendingState: # Radar and noise states precip_cascades: np.ndarray | None = None - precip_nowcast_cascades: np.ndarray | None = None precip_noise_input: np.ndarray | None = None precip_noise_cascades: np.ndarray | None = None precip_mean_noise: np.ndarray | None = None @@ -361,6 +360,8 @@ class StepsBlendingState: # Extrapolation states mean_extrapolation: np.ndarray | None = None std_extrapolation: np.ndarray | None = None + mean_nowcast_timestep: np.ndarray | None = None + std_nowcast_timestep: np.ndarray | None = None rho_extrap_cascade_prev: np.ndarray | None = None rho_extrap_cascade: np.ndarray | None = None precip_cascades_prev_subtimestep: np.ndarray | None = None @@ -1111,6 +1112,16 @@ def transform_to_lagrangian(precip, i): ) precip_forecast_decomp.append(precip_forecast) + # Rearrange the cascaded into a four-dimensional array of shape + # (n_cascade_levels,ar_order+1,m,n) for the autoregressive model + self.__state.precip_cascades = nowcast_utils.stack_cascades( + precip_forecast_decomp, self.__config.n_cascade_levels + ) + + precip_forecast_decomp = precip_forecast_decomp[-1] + self.__state.mean_extrapolation = np.array(precip_forecast_decomp["means"]) + self.__state.std_extrapolation = np.array(precip_forecast_decomp["stds"]) + # Decompose precomputed nowcasts and rearange them again into the required components # TODO: check if this also works if an ensemble nowcast is provided! if self.__precip_nowcast is not None: @@ -1134,60 +1145,63 @@ def transform_to_lagrangian(precip, i): # self.__state.precip_nowcast_cascades = nowcast_utils.stack_cascades( # precip_nowcast_decomp, self.__config.n_cascade_levels # ) - if self.__precip_nowcast.shape[0] == 1: - precip_nowcast_decomp_timestep = [ - self.__params.decomposition_method( - self.__precip_nowcast[0, t, :, :], - bp_filter=self.__params.bandpass_filter, - n_levels=self.__config.n_cascade_levels, - mask=self.__params.mask_threshold, + def decompose_member(member_field, self_obj): + """Loop over timesteps for a single ensemble member.""" + results_levels = [] + results_means = [] + results_stds = [] + for t in range(member_field.shape[0]): # loop over timesteps + res = self_obj.__params.decomposition_method( + field=member_field[t, :, :], + bp_filter=self_obj.__params.bandpass_filter, + n_levels=self_obj.__config.n_cascade_levels, + mask=self_obj.__params.mask_threshold, method="fft", - fft_method=self.__params.fft, - output_domain=self.__config.domain, + fft_method=self_obj.__params.fft, + output_domain=self_obj.__config.domain, compute_stats=True, normalize=True, compact_output=True, ) + results_levels.append(res["cascade_levels"]) + results_means.append(res["means"]) + results_stds.append(res["stds"]) + results = { + "cascade_levels": results_levels, + "means": results_means, + "stds": results_stds, + } + return results + + if self.__precip_nowcast.shape[0] == 1: + precip_nowcast_decomp = [ + decompose_member(self.__precip_nowcast[0], self) ] else: with ThreadPool(self.__config.num_workers) as pool: - precip_nowcast_decomp_timestep = pool.map( - partial( - self.__params.decomposition_method, - bp_filter=self.__params.bandpass_filter, - fft_method=self.__params.fft, - output_domain=self.__config.domain, - normalize=True, - compute_stats=True, - compact_output=True, - ), - list(self.__precip_nowcast[:, t, :, :]), + precip_nowcast_decomp = pool.map( + partial(decompose_member, self_obj=self), + list( + self.__precip_nowcast + ), # each entry = one ensemble member (shape [n_timesteps, x, y]) ) - precip_nowcast_cascades = np.array( - [decomp["cascade_levels"] for decomp in precip_nowcast_decomp_timestep] - ) + self.__state.precip_nowcast_cascades = np.array( + [decomp["cascade_levels"] for decomp in precip_nowcast_decomp] + ).swapaxes(1, 2) - self.__state.precip_nowcast_decomp = nowcast_utils.stack_cascades( - precip_nowcast_cascades, self.__config.n_cascade_levels + print( + "decomposed nowcast shape: ", self.__state.precip_nowcast_cascades.shape ) - # TODO: check if this is also necessary for ensemble nowcasts somwhere, otherwise remove it - precip_nowcast_cascades = precip_nowcast_cascades[-1] self.__state.mean_nowcast = np.array( - np.array(precip_nowcast_cascades["means"]) + [decomp["means"] for decomp in precip_nowcast_decomp] ) - self.__state.std_nowcast = np.array(precip_nowcast_cascades["stds"]) - - # Rearrange the cascaded into a four-dimensional array of shape - # (n_cascade_levels,ar_order+1,m,n) for the autoregressive model - self.__state.precip_cascades = nowcast_utils.stack_cascades( - precip_forecast_decomp, self.__config.n_cascade_levels - ) - - precip_forecast_decomp = precip_forecast_decomp[-1] - self.__state.mean_extrapolation = np.array(precip_forecast_decomp["means"]) - self.__state.std_extrapolation = np.array(precip_forecast_decomp["stds"]) + print("decomposed means shape: ", self.__state.mean_nowcast.shape) + self.__state.std_nowcast = np.array( + [decomp["stds"] for decomp in precip_nowcast_decomp] + ) + print("decomposed stds shape: ", self.__state.std_nowcast.shape) # If necessary, recompose (NWP) model forecasts self.__state.precip_models_cascades = None @@ -1422,7 +1436,8 @@ def __estimate_ar_parameters_radar(self): predefined climatological values. Adjust coefficients if necessary and estimate AR model parameters. """ - + print("mask_threshold is:", self.__params.mask_threshold) + print("cascade shape is:", self.__state.precip_cascades.shape) # If there are values in the radar fields, compute the auto-correlations GAMMA = np.empty((self.__config.n_cascade_levels, self.__config.ar_order)) if not self.__params.zero_precip_radar: @@ -2207,9 +2222,11 @@ def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): # ) for i in range(self.__config.n_cascade_levels): # Use a deterministic Externally computed nowcasting model - worker_state.precip_cascades[j][i] = self.__state.precip_nowcast_decomp[ - j - ][i][t] + worker_state.precip_cascades[j][i] = ( + self.__state.precip_nowcast_cascades[j][i][t] + ) + # worker_state.mean_nowcast_timestep = self.__state.mean_nowcast[:][t] + # worker_state.std_nowcast_timestep = self.__state.std_nowcast[:][t] # Follow the 'standard' STEPS blending approach as described in :cite:`Imhoff2023` elif self.__config.nowcasting_method == "steps": @@ -2676,16 +2693,17 @@ def __blend_cascades(self, t_sub, j, worker_state): ), axis=0, ) # [(extr_field, n_model_fields), n_cascade_levels, ...] + means_stacked = np.concatenate( ( - worker_state.mean_nowcast[None, :], + worker_state.mean_nowcast_timestep, worker_state.mean_models_timestep, ), axis=0, ) sigmas_stacked = np.concatenate( ( - worker_state.std_nowcast[None, :], + worker_state.std_nowcast_timestep, worker_state.std_models_timestep, ), axis=0, @@ -2732,17 +2750,21 @@ def __blend_cascades(self, t_sub, j, worker_state): ), axis=0, ) # [(extr_field, n_model_fields), n_cascade_levels, ...] + print("mean_nowcast timestep shape:", worker_state.mean_nowcast) + print("mean models timestep shape:", worker_state.mean_models_timestep) means_stacked = np.concatenate( ( - worker_state.mean_nowcast[None, :], - worker_state.mean_models_timestep, + # worker_state.mean_nowcast_timestep, + worker_state.mean_extrapolation[None, :], + worker_state.mean_models_timestep[None, j], ), axis=0, ) sigmas_stacked = np.concatenate( ( - worker_state.std_nowcast[None, :], - worker_state.std_models_timestep, + # worker_state.std_nowcast_timestep, + worker_state.std_extrapolation[None, :], + worker_state.std_models_timestep[None, j], ), axis=0, ) @@ -2805,17 +2827,23 @@ def __blend_cascades(self, t_sub, j, worker_state): # Create weights_with_noise to ensure there is always a 3D weights field, even # if self.__config.nowcasting_method is "external_nowcast" and n_ens_members is 1. worker_state.weights_with_noise = worker_state.weights.copy() + print("weights_with_noise:", worker_state.weights_with_noise) worker_state.weights_model_only_with_noise = ( worker_state.weights_model_only.copy() ) if ( - self.__config.nowcasting_method == "external_nowcast" - and self.__config.n_ens_members == 1 + self.__config.nowcasting_method + == "external_nowcast" + # and self.__config.n_ens_members == 1 ): # First determine the weights without noise worker_state.weights = worker_state.weights[:-1, :] / np.sum( worker_state.weights[:-1, :], axis=0 ) + print("weights_without_noise:", worker_state.weights) + print("means_stacked:", means_stacked) + print("std_stacked:", sigmas_stacked) + worker_state.weights_model_only = worker_state.weights_model_only[ :-1, : ] / np.sum(worker_state.weights_model_only[:-1, :], axis=0) From 9c9154fd10fdf0f6cf912ab23f154b7766522d5b Mon Sep 17 00:00:00 2001 From: Joep1999 Date: Wed, 17 Sep 2025 14:03:27 +0200 Subject: [PATCH 04/18] removed prints, mostly working --- pysteps/blending/steps.py | 78 ++++----------------------------------- 1 file changed, 8 insertions(+), 70 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 95faeb8dc..98e2ec3ef 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -1125,31 +1125,10 @@ def transform_to_lagrangian(precip, i): # Decompose precomputed nowcasts and rearange them again into the required components # TODO: check if this also works if an ensemble nowcast is provided! if self.__precip_nowcast is not None: - # precip_nowcast_decomp = [] - # for i in range(self.__precip_nowcast.shape[0]): - # cascades = self.__params.decomposition_method( - # self.__precip_nowcast[i], - # self.__params.bandpass_filter, - # n_levels=self.__config.n_cascade_levels, - # mask=self.__params.mask_threshold, - # method="fft", - # fft_method=self.__params.fft, - # output_domain=self.__config.domain, - # compute_stats=True, - # normalize=True, - # compact_output=True, - # ) - - # precip_nowcast_decomp.append(cascades) - - # self.__state.precip_nowcast_cascades = nowcast_utils.stack_cascades( - # precip_nowcast_decomp, self.__config.n_cascade_levels - # ) + def decompose_member(member_field, self_obj): """Loop over timesteps for a single ensemble member.""" - results_levels = [] - results_means = [] - results_stds = [] + results = [] for t in range(member_field.shape[0]): # loop over timesteps res = self_obj.__params.decomposition_method( field=member_field[t, :, :], @@ -1163,19 +1142,13 @@ def decompose_member(member_field, self_obj): normalize=True, compact_output=True, ) - results_levels.append(res["cascade_levels"]) - results_means.append(res["means"]) - results_stds.append(res["stds"]) - results = { - "cascade_levels": results_levels, - "means": results_means, - "stds": results_stds, - } + results.append(res["cascade_levels"]) + return results if self.__precip_nowcast.shape[0] == 1: precip_nowcast_decomp = [ - decompose_member(self.__precip_nowcast[0], self) + decompose_member(self.__precip_nowcast[0], self_obj=self) ] else: with ThreadPool(self.__config.num_workers) as pool: @@ -1187,22 +1160,9 @@ def decompose_member(member_field, self_obj): ) self.__state.precip_nowcast_cascades = np.array( - [decomp["cascade_levels"] for decomp in precip_nowcast_decomp] + [precip_nowcast_decomp] ).swapaxes(1, 2) - print( - "decomposed nowcast shape: ", self.__state.precip_nowcast_cascades.shape - ) - - self.__state.mean_nowcast = np.array( - [decomp["means"] for decomp in precip_nowcast_decomp] - ) - print("decomposed means shape: ", self.__state.mean_nowcast.shape) - self.__state.std_nowcast = np.array( - [decomp["stds"] for decomp in precip_nowcast_decomp] - ) - print("decomposed stds shape: ", self.__state.std_nowcast.shape) - # If necessary, recompose (NWP) model forecasts self.__state.precip_models_cascades = None @@ -1436,8 +1396,6 @@ def __estimate_ar_parameters_radar(self): predefined climatological values. Adjust coefficients if necessary and estimate AR model parameters. """ - print("mask_threshold is:", self.__params.mask_threshold) - print("cascade shape is:", self.__state.precip_cascades.shape) # If there are values in the radar fields, compute the auto-correlations GAMMA = np.empty((self.__config.n_cascade_levels, self.__config.ar_order)) if not self.__params.zero_precip_radar: @@ -2216,17 +2174,11 @@ def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): # If nowcast method seleced is external_nowcast, n_ens members has to be 1 if self.__config.nowcasting_method == "external_nowcast": print("Using nowcasting method:", self.__config.nowcasting_method) - # if self.__config.n_ens_members != 1: - # raise ValueError( - # "Currently only 1 ensemble member supported when using (deterministic) external_nowcast as input." - # ) for i in range(self.__config.n_cascade_levels): # Use a deterministic Externally computed nowcasting model worker_state.precip_cascades[j][i] = ( self.__state.precip_nowcast_cascades[j][i][t] ) - # worker_state.mean_nowcast_timestep = self.__state.mean_nowcast[:][t] - # worker_state.std_nowcast_timestep = self.__state.std_nowcast[:][t] # Follow the 'standard' STEPS blending approach as described in :cite:`Imhoff2023` elif self.__config.nowcasting_method == "steps": @@ -2294,11 +2246,7 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( worker_state.noise_extrapolated_decomp = [] worker_state.precip_extrapolated_probability_matching = [] - if ( - self.__config.nowcasting_method - == "external_nowcast" - # and self.__config.n_ens_members == 1 - ): + if self.__config.nowcasting_method == "external_nowcast": for i in range(self.__config.n_cascade_levels): precip_extrapolated_decomp = worker_state.precip_cascades[j][i][-1, :] worker_state.time_prev_timestep[j] = t + 1 @@ -2750,8 +2698,6 @@ def __blend_cascades(self, t_sub, j, worker_state): ), axis=0, ) # [(extr_field, n_model_fields), n_cascade_levels, ...] - print("mean_nowcast timestep shape:", worker_state.mean_nowcast) - print("mean models timestep shape:", worker_state.mean_models_timestep) means_stacked = np.concatenate( ( # worker_state.mean_nowcast_timestep, @@ -2827,22 +2773,14 @@ def __blend_cascades(self, t_sub, j, worker_state): # Create weights_with_noise to ensure there is always a 3D weights field, even # if self.__config.nowcasting_method is "external_nowcast" and n_ens_members is 1. worker_state.weights_with_noise = worker_state.weights.copy() - print("weights_with_noise:", worker_state.weights_with_noise) worker_state.weights_model_only_with_noise = ( worker_state.weights_model_only.copy() ) - if ( - self.__config.nowcasting_method - == "external_nowcast" - # and self.__config.n_ens_members == 1 - ): + if self.__config.nowcasting_method == "external_nowcast": # First determine the weights without noise worker_state.weights = worker_state.weights[:-1, :] / np.sum( worker_state.weights[:-1, :], axis=0 ) - print("weights_without_noise:", worker_state.weights) - print("means_stacked:", means_stacked) - print("std_stacked:", sigmas_stacked) worker_state.weights_model_only = worker_state.weights_model_only[ :-1, : From 06c648304725c223cf93482c1c0e16d8cb181613 Mon Sep 17 00:00:00 2001 From: Ruben Imhoff Date: Wed, 17 Sep 2025 17:29:10 +0200 Subject: [PATCH 05/18] fix: make sure means and stds of external nowcast cascade are used --- pysteps/blending/steps.py | 580 +++++++++++++------------------------- 1 file changed, 198 insertions(+), 382 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 98e2ec3ef..26869b12b 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -352,6 +352,7 @@ class StepsBlendingParams: class StepsBlendingState: # Radar and noise states precip_cascades: np.ndarray | None = None + precip_nowcast_cascades: np.ndarray | None = None precip_noise_input: np.ndarray | None = None precip_noise_cascades: np.ndarray | None = None precip_mean_noise: np.ndarray | None = None @@ -571,9 +572,7 @@ def compute_forecast(self): self.__prepare_forecast_loop() self.__initialize_noise_cascades() if self.__config.measure_time: - self.__init_time = self.__measure_time( - "initialization", self.__start_time_init - ) + self.__init_time = self.__measure_time("initialization", self.__start_time_init) self.__blended_nowcast_main_loop() # Stack and return the forecast output @@ -608,16 +607,10 @@ def __blended_nowcast_main_loop(self): starttime_mainloop = time.time() self.__state.extrapolation_kwargs["return_displacement"] = True - self.__state.precip_cascades_prev_subtimestep = deepcopy( - self.__state.precip_cascades - ) - self.__state.cascade_noise_prev_subtimestep = deepcopy( - self.__state.precip_noise_cascades - ) + self.__state.precip_cascades_prev_subtimestep = deepcopy(self.__state.precip_cascades) + self.__state.cascade_noise_prev_subtimestep = deepcopy(self.__state.precip_noise_cascades) - self.__state.time_prev_timestep = [ - 0.0 for j in range(self.__config.n_ens_members) - ] + self.__state.time_prev_timestep = [0.0 for j in range(self.__config.n_ens_members)] self.__state.leadtime_since_start_forecast = [ 0.0 for j in range(self.__config.n_ens_members) ] @@ -650,13 +643,11 @@ def worker(j): if t_sub > 0: self.__blend_cascades(t_sub, j, worker_state) self.__recompose_cascade_to_rainfall_field(j, worker_state) - final_blended_forecast_single_member = ( - self.__post_process_output( - j, - t_sub, - final_blended_forecast_single_member, - worker_state, - ) + final_blended_forecast_single_member = self.__post_process_output( + j, + t_sub, + final_blended_forecast_single_member, + worker_state, ) final_blended_forecast_all_members_one_timestep[j] = ( final_blended_forecast_single_member @@ -684,9 +675,7 @@ def worker(j): print("done.") if self.__config.callback is not None: - precip_forecast_final = np.stack( - final_blended_forecast_all_members_one_timestep - ) + precip_forecast_final = np.stack(final_blended_forecast_all_members_one_timestep) if precip_forecast_final.shape[1] > 0: self.__config.callback(precip_forecast_final.squeeze()) @@ -724,19 +713,14 @@ def __check_inputs(self): raise KeyError( "if precip_nowcast is not None, nowcasting_method should be set to 'external_nowcast' " ) - if ( - self.__config.nowcasting_method == "external_nowcast" - and self.__precip_nowcast is None - ): + if self.__config.nowcasting_method == "external_nowcast" and self.__precip_nowcast is None: raise KeyError( "if nowcasting_method is set to 'external_nowcast', an external precip_nowcast should be provided as variable." ) # Check dimensions of velocity if self.__velocity.ndim != 3: - raise ValueError( - "velocity must be a three-dimensional array of shape (2, m, n)" - ) + raise ValueError("velocity must be a three-dimensional array of shape (2, m, n)") if self.__velocity_models.ndim != 5: raise ValueError( "velocity_models must be a five-dimensional array of shape (n_models, timestep, 2, m, n)" @@ -761,17 +745,13 @@ def __check_inputs(self): if isinstance(self.__timesteps, list): self.__params.time_steps_is_list = True if not sorted(self.__timesteps) == self.__timesteps: - raise ValueError( - "timesteps is not in ascending order", self.__timesteps - ) + raise ValueError("timesteps is not in ascending order", self.__timesteps) if self.__precip_models.shape[1] != math.ceil(self.__timesteps[-1]) + 1: raise ValueError( "precip_models does not contain sufficient lead times for this forecast" ) self.__params.original_timesteps = [0] + list(self.__timesteps) - self.__timesteps = nowcast_utils.binned_timesteps( - self.__params.original_timesteps - ) + self.__timesteps = nowcast_utils.binned_timesteps(self.__params.original_timesteps) else: self.__params.time_steps_is_list = False if self.__precip_models.shape[1] != self.__timesteps + 1: @@ -804,9 +784,7 @@ def __check_inputs(self): if self.__config.extrapolation_kwargs is None: self.__state.extrapolation_kwargs = dict() else: - self.__state.extrapolation_kwargs = deepcopy( - self.__config.extrapolation_kwargs - ) + self.__state.extrapolation_kwargs = deepcopy(self.__config.extrapolation_kwargs) if self.__config.filter_kwargs is None: self.__params.filter_kwargs = dict() @@ -827,13 +805,9 @@ def __check_inputs(self): if self.__config.climatology_kwargs is None: # Make sure clim_kwargs at least contains the number of models - self.__params.climatology_kwargs = dict( - {"n_models": self.__precip_models.shape[0]} - ) + self.__params.climatology_kwargs = dict({"n_models": self.__precip_models.shape[0]}) else: - self.__params.climatology_kwargs = deepcopy( - self.__config.climatology_kwargs - ) + self.__params.climatology_kwargs = deepcopy(self.__config.climatology_kwargs) if self.__config.mask_kwargs is None: self.__params.mask_kwargs = dict() @@ -854,10 +828,7 @@ def __check_inputs(self): if self.__config.conditional and self.__params.precip_threshold is None: raise ValueError("conditional=True but precip_thr is not set") - if ( - self.__config.mask_method is not None - and self.__params.precip_threshold is None - ): + if self.__config.mask_method is not None and self.__params.precip_threshold is None: raise ValueError("mask_method!=None but precip_thr=None") if self.__config.noise_stddev_adj not in ["auto", "fixed", None]: @@ -868,17 +839,13 @@ def __check_inputs(self): if self.__config.kmperpixel is None: if self.__config.velocity_perturbation_method is not None: - raise ValueError( - "velocity_perturbation_method is set but kmperpixel=None" - ) + raise ValueError("velocity_perturbation_method is set but kmperpixel=None") if self.__config.mask_method == "incremental": raise ValueError("mask_method='incremental' but kmperpixel=None") if self.__config.timestep is None: if self.__config.velocity_perturbation_method is not None: - raise ValueError( - "velocity_perturbation_method is set but timestep=None" - ) + raise ValueError("velocity_perturbation_method is set but timestep=None") if self.__config.mask_method == "incremental": raise ValueError("mask_method='incremental' but timestep=None") @@ -893,9 +860,7 @@ def __print_forecast_info(self): print("Inputs") print("------") print(f"forecast issue time: {self.__issuetime.isoformat()}") - print( - f"input dimensions: {self.__precip.shape[1]}x{self.__precip.shape[2]}" - ) + print(f"input dimensions: {self.__precip.shape[1]}x{self.__precip.shape[2]}") if self.__precip_nowcast is not None: print( f"input dimensions pre-computed nowcast: {self.__precip_nowcast.shape[2]}x{self.__precip_nowcast.shape[3]}" @@ -910,9 +875,7 @@ def __print_forecast_info(self): print("-----------------------") print(f"number of (NWP) models: {self.__precip_models.shape[0]}") print(f"blend (NWP) model members: {self.__config.blend_nwp_members}") - print( - f"decompose (NWP) models: {'yes' if self.__precip_models.ndim == 4 else 'no'}" - ) + print(f"decompose (NWP) models: {'yes' if self.__precip_models.ndim == 4 else 'no'}") print("") print("Methods") @@ -922,16 +885,10 @@ def __print_forecast_info(self): print(f"decomposition: {self.__config.decomposition_method}") print(f"nowcasting algorithm: {self.__config.nowcasting_method}") print(f"noise generator: {self.__config.noise_method}") - print( - f"noise adjustment: {'yes' if self.__config.noise_stddev_adj else 'no'}" - ) - print( - f"velocity perturbator: {self.__config.velocity_perturbation_method}" - ) + print(f"noise adjustment: {'yes' if self.__config.noise_stddev_adj else 'no'}") + print(f"velocity perturbator: {self.__config.velocity_perturbation_method}") print(f"blending weights method: {self.__config.weights_method}") - print( - f"conditional statistics: {'yes' if self.__config.conditional else 'no'}" - ) + print(f"conditional statistics: {'yes' if self.__config.conditional else 'no'}") print(f"precip. mask method: {self.__config.mask_method}") print(f"probability matching: {self.__config.probmatching_method}") print(f"FFT method: {self.__config.fft_method}") @@ -1071,29 +1028,24 @@ def transform_to_lagrangian(precip, i): for i in range(self.__config.ar_order): res.append(dask.delayed(transform_to_lagrangian)(self.__precip, i)) num_workers_ = ( - len(res) - if self.__config.num_workers > len(res) - else self.__config.num_workers + len(res) if self.__config.num_workers > len(res) else self.__config.num_workers ) self.__precip = np.stack( - list(dask.compute(*res, num_workers=num_workers_)) - + [self.__precip[-1, :, :]] + list(dask.compute(*res, num_workers=num_workers_)) + [self.__precip[-1, :, :]] ) # Replace non-finite values with the minimum value for each field self.__precip = self.__precip.copy() for i in range(self.__precip.shape[0]): - self.__precip[i, ~np.isfinite(self.__precip[i, :])] = np.nanmin( - self.__precip[i, :] - ) + self.__precip[i, ~np.isfinite(self.__precip[i, :])] = np.nanmin(self.__precip[i, :]) if self.__precip_nowcast is not None: self.__precip_nowcast = self.__precip_nowcast.copy() # TODO: check this function for ensemble nowcasts -> now made it always work with 4 dimentions for m in range(self.__precip_nowcast.shape[0]): for t in range(self.__precip_nowcast.shape[1]): - self.__precip_nowcast[ - m, t, ~np.isfinite(self.__precip_nowcast[m, t, :, :]) - ] = np.nanmin(self.__precip_nowcast[m, t, :, :]) + self.__precip_nowcast[m, t, ~np.isfinite(self.__precip_nowcast[m, t, :, :])] = ( + np.nanmin(self.__precip_nowcast[m, t, :, :]) + ) # Perform the cascade decomposition for the input precip fields and, # if necessary, for the (NWP) model fields @@ -1123,45 +1075,23 @@ def transform_to_lagrangian(precip, i): self.__state.std_extrapolation = np.array(precip_forecast_decomp["stds"]) # Decompose precomputed nowcasts and rearange them again into the required components - # TODO: check if this also works if an ensemble nowcast is provided! if self.__precip_nowcast is not None: - - def decompose_member(member_field, self_obj): - """Loop over timesteps for a single ensemble member.""" - results = [] - for t in range(member_field.shape[0]): # loop over timesteps - res = self_obj.__params.decomposition_method( - field=member_field[t, :, :], - bp_filter=self_obj.__params.bandpass_filter, - n_levels=self_obj.__config.n_cascade_levels, - mask=self_obj.__params.mask_threshold, - method="fft", - fft_method=self_obj.__params.fft, - output_domain=self_obj.__config.domain, - compute_stats=True, - normalize=True, - compact_output=True, - ) - results.append(res["cascade_levels"]) - - return results - if self.__precip_nowcast.shape[0] == 1: - precip_nowcast_decomp = [ - decompose_member(self.__precip_nowcast[0], self_obj=self) - ] + precip_nowcast_decomp, precip_nowcast_means, precip_nowcast_stds = ( + self.__decompose_member(self.__precip_nowcast[0]) + ) else: with ThreadPool(self.__config.num_workers) as pool: - precip_nowcast_decomp = pool.map( - partial(decompose_member, self_obj=self), + precip_nowcast_decomp, precip_nowcast_means, precip_nowcast_stds = pool.map( + partial(self.__decompose_member), list( self.__precip_nowcast ), # each entry = one ensemble member (shape [n_timesteps, x, y]) ) - self.__state.precip_nowcast_cascades = np.array( - [precip_nowcast_decomp] - ).swapaxes(1, 2) + self.__state.precip_nowcast_cascades = np.array([precip_nowcast_decomp]).swapaxes(1, 2) + self.__state.mean_nowcast_timestep = np.array([precip_nowcast_means]).swapaxes(1, 2) + self.__state.std_nowcast_timestep = np.array([precip_nowcast_stds]).swapaxes(1, 2) # If necessary, recompose (NWP) model forecasts self.__state.precip_models_cascades = None @@ -1200,15 +1130,37 @@ def decompose_member(member_field, self_obj): self.__params.noise_kwargs["win_fun"], ) + def __decompose_member(self, member_field): + """Loop over timesteps for a single ensemble member.""" + results_decomp = [] + means = [] + stds = [] + for t in range(member_field.shape[0]): # loop over timesteps + res = self.__params.decomposition_method( + field=member_field[t, :, :], + bp_filter=self.__params.bandpass_filter, + n_levels=self.__config.n_cascade_levels, + mask=self.__params.mask_threshold, + method="fft", + fft_method=self.__params.fft, + output_domain=self.__config.domain, + compute_stats=True, + normalize=True, + compact_output=True, + ) + results_decomp.append(res["cascade_levels"]) + means.append(res["means"]) + stds.append(res["stds"]) + + return results_decomp, means, stds + def __zero_precipitation_forecast(self): """ Generate a zero-precipitation forecast (filled with the minimum precip value) when no precipitation above the threshold is detected. The forecast is optionally returned or passed to a callback. """ - print( - "No precipitation above the threshold found in both the radar and NWP fields" - ) + print("No precipitation above the threshold found in both the radar and NWP fields") print("The resulting forecast will contain only zeros") # Create the output list precip_forecast = [[] for j in range(self.__config.n_ens_members)] @@ -1239,10 +1191,7 @@ def __zero_precipitation_forecast(self): if self.__config.return_output: precip_forecast_all_members_all_times = np.stack( - [ - np.stack(precip_forecast[j]) - for j in range(self.__config.n_ens_members) - ] + [np.stack(precip_forecast[j]) for j in range(self.__config.n_ens_members)] ) if self.__config.measure_time: @@ -1266,9 +1215,7 @@ def __prepare_nowcast_for_zero_radar(self): # forecast. I.e., velocity (radar) equals velocity_model at the first time # step. # Use the velocity from velocity_models at time step 0 - self.__velocity = self.__velocity_models[:, 0, :, :, :].astype( - np.float64, copy=False - ) + self.__velocity = self.__velocity_models[:, 0, :, :, :].astype(np.float64, copy=False) # Take the average over the first axis, which corresponds to n_models # (hence, the model average) self.__velocity = np.mean(self.__velocity, axis=0) @@ -1294,9 +1241,7 @@ def __prepare_nowcast_for_zero_radar(self): max_rain_pixels = rain_pixels max_rain_pixels_j = j max_rain_pixels_t = t - self.__state.precip_noise_input = self.__precip_models[max_rain_pixels_j][ - max_rain_pixels_t - ] + self.__state.precip_noise_input = self.__precip_models[max_rain_pixels_j][max_rain_pixels_t] self.__state.precip_noise_input = self.__state.precip_noise_input.astype( np.float64, copy=False ) @@ -1304,12 +1249,10 @@ def __prepare_nowcast_for_zero_radar(self): # If zero_precip_radar, make sure that precip_cascade does not contain # only nans or infs. If so, fill it with the zero value. if self.__state.precip_models_cascades is not None: - self.__state.precip_cascades[~np.isfinite(self.__state.precip_cascades)] = ( - np.nanmin( - self.__state.precip_models_cascades[ - max_rain_pixels_j, max_rain_pixels_t - ]["cascade_levels"] - ) + self.__state.precip_cascades[~np.isfinite(self.__state.precip_cascades)] = np.nanmin( + self.__state.precip_models_cascades[max_rain_pixels_j, max_rain_pixels_t][ + "cascade_levels" + ] ) else: precip_models_cascade_timestep = self.__params.decomposition_method( @@ -1321,15 +1264,13 @@ def __prepare_nowcast_for_zero_radar(self): compute_stats=True, compact_output=True, )["cascade_levels"] - self.__state.precip_cascades[~np.isfinite(self.__state.precip_cascades)] = ( - np.nanmin(precip_models_cascade_timestep) + self.__state.precip_cascades[~np.isfinite(self.__state.precip_cascades)] = np.nanmin( + precip_models_cascade_timestep ) # Make sure precip_noise_input is three-dimensional if len(self.__state.precip_noise_input.shape) != 3: - self.__state.precip_noise_input = self.__state.precip_noise_input[ - np.newaxis, :, : - ] + self.__state.precip_noise_input = self.__state.precip_noise_input[np.newaxis, :, :] def __initialize_noise(self): """ @@ -1338,9 +1279,7 @@ def __initialize_noise(self): """ if self.__config.noise_method is not None: # get methods for perturbations - init_noise, self.__params.noise_generator = noise.get_method( - self.__config.noise_method - ) + init_noise, self.__params.noise_generator = noise.get_method(self.__config.noise_method) # initialize the perturbation generator for the precipitation field self.__params.perturbation_generator = init_noise( @@ -1452,15 +1391,11 @@ def __estimate_ar_parameters_radar(self): # adjust the lag-2 correlation coefficient to ensure that the AR(p) # process is stationary for i in range(self.__config.n_cascade_levels): - GAMMA[i, 1] = autoregression.adjust_lag2_corrcoef2( - GAMMA[i, 0], GAMMA[i, 1] - ) + GAMMA[i, 1] = autoregression.adjust_lag2_corrcoef2(GAMMA[i, 0], GAMMA[i, 1]) # estimate the parameters of the AR(p) model from the auto-correlation # coefficients - self.__params.PHI = np.empty( - (self.__config.n_cascade_levels, self.__config.ar_order + 1) - ) + self.__params.PHI = np.empty((self.__config.n_cascade_levels, self.__config.ar_order + 1)) for i in range(self.__config.n_cascade_levels): self.__params.PHI[i, :] = autoregression.estimate_ar_params_yw(GAMMA[i, :]) @@ -1551,16 +1486,12 @@ def __prepare_forecast_loop(self): self.__state.previous_displacement_prob_matching = np.stack( [None for j in range(self.__config.n_ens_members)] ) - self.__state.final_blended_forecast = [ - [] for j in range(self.__config.n_ens_members) - ] + self.__state.final_blended_forecast = [[] for j in range(self.__config.n_ens_members)] if self.__config.mask_method == "incremental": # get mask parameters self.__params.mask_rim = self.__params.mask_kwargs.get("mask_rim", 10) - self.__params.max_mask_rim = self.__params.mask_kwargs.get( - "max_mask_rim", 10 - ) + self.__params.max_mask_rim = self.__params.mask_kwargs.get("max_mask_rim", 10) mask_f = self.__params.mask_kwargs.get("mask_f", 1.0) # initialize the structuring element struct = generate_binary_structure(2, 1) @@ -1589,12 +1520,8 @@ def __prepare_forecast_loop(self): # initizalize the current and previous extrapolation forecast scale for the nowcasting component # phi1 / (1 - phi2), see BPS2004 - self.__state.rho_extrap_cascade_prev = np.repeat( - 1.0, self.__params.PHI.shape[0] - ) - self.__state.rho_extrap_cascade = self.__params.PHI[:, 0] / ( - 1.0 - self.__params.PHI[:, 1] - ) + self.__state.rho_extrap_cascade_prev = np.repeat(1.0, self.__params.PHI.shape[0]) + self.__state.rho_extrap_cascade = self.__params.PHI[:, 0] / (1.0 - self.__params.PHI[:, 1]) def __initialize_noise_cascades(self): """ @@ -1602,9 +1529,7 @@ def __initialize_noise_cascades(self): We also need to return the mean and standard deviations of the noise for the recombination of the noise before advecting it. """ - self.__state.precip_noise_cascades = np.zeros( - self.__state.precip_cascades.shape - ) + self.__state.precip_noise_cascades = np.zeros(self.__state.precip_cascades.shape) self.__state.precip_mean_noise = np.zeros( (self.__config.n_ens_members, self.__config.n_cascade_levels) ) @@ -1725,17 +1650,13 @@ def __decompose_nwp_if_needed_and_fill_nans_in_nwp(self, t): self.__state.precip_models_cascades_timestep[ ~np.isfinite(self.__state.precip_models_cascades_timestep) ] = min_cascade - self.__state.precip_models_timestep[ - ~np.isfinite(self.__state.precip_models_timestep) - ] = min_precip + self.__state.precip_models_timestep[~np.isfinite(self.__state.precip_models_timestep)] = ( + min_precip + ) # Also set any nans or infs in the mean and sigma of the cascade to # respectively 0.0 and 1.0 - self.__state.mean_models_timestep[ - ~np.isfinite(self.__state.mean_models_timestep) - ] = 0.0 - self.__state.std_models_timestep[ - ~np.isfinite(self.__state.std_models_timestep) - ] = 0.0 + self.__state.mean_models_timestep[~np.isfinite(self.__state.mean_models_timestep)] = 0.0 + self.__state.std_models_timestep[~np.isfinite(self.__state.std_models_timestep)] = 0.0 def __find_nowcast_NWP_combination(self, t): """ @@ -1747,9 +1668,9 @@ def __find_nowcast_NWP_combination(self, t): np.float64, copy=False ) - self.__state.velocity_models_timestep = self.__velocity_models[ - :, t, :, :, : - ].astype(np.float64, copy=False) + self.__state.velocity_models_timestep = self.__velocity_models[:, t, :, :, :].astype( + np.float64, copy=False + ) # Make sure the number of model members is not larger than or equal to n_ens_members n_model_members = self.__state.precip_models_cascades_timestep.shape[0] if n_model_members > self.__config.n_ens_members: @@ -1777,9 +1698,7 @@ def __find_nowcast_NWP_combination(self, t): # for indexing the right climatological skill file when pysteps calculates # the blended forecast in parallel. if n_model_members > 1: - self.__state.mapping_list_NWP_member_to_ensemble_member = np.arange( - n_model_members - ) + self.__state.mapping_list_NWP_member_to_ensemble_member = np.arange(n_model_members) else: self.__state.mapping_list_NWP_member_to_ensemble_member = [0] @@ -1825,8 +1744,8 @@ def __find_nowcast_NWP_combination(self, t): n_ens_members_max, axis=0, ) - self.__state.mean_nowcast = np.repeat( - self.__state.mean_nowcast, n_ens_members_max, axis=0 + self.__state.mean_nowcast_timestep = np.repeat( + self.__state.mean_nowcast_timestep, n_ens_members_max, axis=0 ) self.__state.std_nowcast_timestep = np.repeat( self.__state.std_nowcast, n_ens_members_max, axis=0 @@ -1880,12 +1799,10 @@ def __find_nowcast_NWP_combination(self, t): repeats, axis=0, ) - self.__state.mean_nowcast = np.repeat( - self.__state.mean_nowcast, repeats, axis=0 - ) - self.__state.std_nowcast = np.repeat( - self.__state.std_nowcast, repeats, axis=0 + self.__state.mean_nowcast_timestep = np.repeat( + self.__state.mean_nowcast_timestep, repeats, axis=0 ) + self.__state.std_nowcast = np.repeat(self.__state.std_nowcast, repeats, axis=0) self.__state.velocity_nowcast = np.repeat( self.__state.velocity_nowcast, repeats, axis=0 ) @@ -1903,9 +1820,7 @@ def __find_nowcast_NWP_combination(self, t): # for indexing the right climatological skill file when pysteps calculates # the blended forecast in parallel. if n_model_members > 1: - self.__state.mapping_list_NWP_member_to_ensemble_member = np.arange( - n_model_members - ) + self.__state.mapping_list_NWP_member_to_ensemble_member = np.arange(n_model_members) else: self.__state.mapping_list_NWP_member_to_ensemble_member = [0] @@ -1983,34 +1898,26 @@ def __determine_skill_for_current_timestep(self, t): if t == 0: # Calculate the initial skill of the (NWP) model forecasts at t=0. self.__params.rho_nwp_models = [] - for model_index in range( - self.__state.precip_models_cascades_timestep.shape[0] - ): + for model_index in range(self.__state.precip_models_cascades_timestep.shape[0]): rho_value = blending.skill_scores.spatial_correlation( obs=self.__state.precip_cascades[0, :, -1, :, :].copy(), - mod=self.__state.precip_models_cascades_timestep[ - model_index, :, :, : - ].copy(), + mod=self.__state.precip_models_cascades_timestep[model_index, :, :, :].copy(), domain_mask=self.__params.domain_mask, ) self.__params.rho_nwp_models.append(rho_value) self.__params.rho_nwp_models = np.stack(self.__params.rho_nwp_models) # Ensure that the model skill decreases with increasing scale level. - for model_index in range( - self.__state.precip_models_cascades_timestep.shape[0] - ): - for i in range( - 1, self.__state.precip_models_cascades_timestep.shape[1] - ): + for model_index in range(self.__state.precip_models_cascades_timestep.shape[0]): + for i in range(1, self.__state.precip_models_cascades_timestep.shape[1]): if ( self.__params.rho_nwp_models[model_index, i] > self.__params.rho_nwp_models[model_index, i - 1] ): # Set it equal to the previous scale level - self.__params.rho_nwp_models[model_index, i] = ( - self.__params.rho_nwp_models[model_index, i - 1] - ) + self.__params.rho_nwp_models[model_index, i] = self.__params.rho_nwp_models[ + model_index, i - 1 + ] # Save this in the climatological skill file blending.clim.save_skill( @@ -2079,9 +1986,7 @@ def __determine_weights_per_component(self, worker_state): # selected, weights will be overwritten with those weights prior to # blending step. # weight = [(extr_field, n_model_fields, noise), n_cascade_levels, ...] - worker_state.weights = calculate_weights_bps( - worker_state.rho_final_blended_forecast - ) + worker_state.weights = calculate_weights_bps(worker_state.rho_final_blended_forecast) # The model only weights if self.__config.weights_method == "bps": @@ -2114,9 +2019,7 @@ def __determine_weights_per_component(self, worker_state): n_model, i, :, : ].flatten() for n_model in range( - worker_state.precip_models_cascades_timestep.shape[ - 0 - ] + worker_state.precip_models_cascades_timestep.shape[0] ) ] ) @@ -2133,8 +2036,7 @@ def __determine_weights_per_component(self, worker_state): ) else: raise ValueError( - "Unknown weights method %s: must be 'bps' or 'spn'" - % self.__config.weights_method + "Unknown weights method %s: must be 'bps' or 'spn'" % self.__config.weights_method ) def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): @@ -2176,9 +2078,7 @@ def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): print("Using nowcasting method:", self.__config.nowcasting_method) for i in range(self.__config.n_cascade_levels): # Use a deterministic Externally computed nowcasting model - worker_state.precip_cascades[j][i] = ( - self.__state.precip_nowcast_cascades[j][i][t] - ) + worker_state.precip_cascades[j][i] = self.__state.precip_nowcast_cascades[j][i][t] # Follow the 'standard' STEPS blending approach as described in :cite:`Imhoff2023` elif self.__config.nowcasting_method == "steps": @@ -2189,10 +2089,8 @@ def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): epsilon_decomposed is not None or self.__config.velocity_perturbation_method is not None ): - worker_state.precip_cascades[j][i] = ( - autoregression.iterate_ar_model( - worker_state.precip_cascades[j][i], self.__params.PHI[i, :] - ) + worker_state.precip_cascades[j][i] = autoregression.iterate_ar_model( + worker_state.precip_cascades[j][i], self.__params.PHI[i, :] ) # Renormalize the cascade worker_state.precip_cascades[j][i][1] /= np.std( @@ -2217,12 +2115,10 @@ def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): epsilon_temp = None # apply AR(p) process to noise cascade level # (Returns zero noise if epsilon_decomposed is None) - worker_state.precip_noise_cascades[j][i] = ( - autoregression.iterate_ar_model( - worker_state.precip_noise_cascades[j][i], - self.__params.PHI[i, :], - eps=epsilon_temp, - ) + worker_state.precip_noise_cascades[j][i] = autoregression.iterate_ar_model( + worker_state.precip_noise_cascades[j][i], + self.__params.PHI[i, :], + eps=epsilon_temp, ) epsilon_decomposed = None @@ -2251,9 +2147,7 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( precip_extrapolated_decomp = worker_state.precip_cascades[j][i][-1, :] worker_state.time_prev_timestep[j] = t + 1 - worker_state.precip_extrapolated_decomp.append( - precip_extrapolated_decomp.copy() - ) + worker_state.precip_extrapolated_decomp.append(precip_extrapolated_decomp.copy()) worker_state.precip_extrapolated_probability_matching.append( precip_extrapolated_decomp.copy() ) @@ -2303,9 +2197,7 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( noise_cascade_subtimestep = np.stack(noise_cascade_subtimestep) t_diff_prev_subtimestep = t_sub - worker_state.time_prev_timestep[j] - worker_state.leadtime_since_start_forecast[ - j - ] += t_diff_prev_subtimestep + worker_state.leadtime_since_start_forecast[j] += t_diff_prev_subtimestep # compute the perturbed motion field - include the NWP # velocities and the weights. Note that we only perturb @@ -2355,12 +2247,10 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( # This is needed to remove the interpolation artifacts. # In addition, the number of extrapolations is greatly reduced # A. Radar Rain - precip_forecast_recomp_subtimestep = ( - blending.utils.recompose_cascade( - combined_cascade=precip_forecast_cascade_subtimestep, - combined_mean=worker_state.mean_extrapolation, - combined_sigma=worker_state.std_extrapolation, - ) + precip_forecast_recomp_subtimestep = blending.utils.recompose_cascade( + combined_cascade=precip_forecast_cascade_subtimestep, + combined_mean=worker_state.mean_extrapolation, + combined_sigma=worker_state.std_extrapolation, ) # Make sure we have values outside the mask if self.__params.zero_precip_radar: @@ -2372,9 +2262,7 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( neginf=self.__params.precip_zerovalue, ) # Put back the mask - precip_forecast_recomp_subtimestep[self.__params.domain_mask] = ( - np.nan - ) + precip_forecast_recomp_subtimestep[self.__params.domain_mask] = np.nan worker_state.extrapolation_kwargs["displacement_prev"] = ( worker_state.previous_displacement[j] ) @@ -2449,17 +2337,13 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( compact_output=True, )["cascade_levels"] for i in range(self.__config.n_cascade_levels): - noise_extrapolated_decomp[i] *= self.__params.noise_std_coeffs[ - i - ] + noise_extrapolated_decomp[i] *= self.__params.noise_std_coeffs[i] # Append the results to the output lists worker_state.precip_extrapolated_decomp.append( precip_extrapolated_decomp.copy() ) - worker_state.noise_extrapolated_decomp.append( - noise_extrapolated_decomp.copy() - ) + worker_state.noise_extrapolated_decomp.append(noise_extrapolated_decomp.copy()) precip_forecast_cascade_subtimestep = None precip_forecast_recomp_subtimestep = None precip_forecast_extrapolated_recomp_subtimestep_temp = None @@ -2480,9 +2364,9 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( ) # Apply the domain mask to the extrapolation component precip_forecast_temp_for_probability_matching = self.__precip.copy() - precip_forecast_temp_for_probability_matching[ - self.__params.domain_mask - ] = np.nan + precip_forecast_temp_for_probability_matching[self.__params.domain_mask] = ( + np.nan + ) ( precip_forecast_extrapolated_probability_matching_temp, @@ -2525,8 +2409,7 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( self.__velocity + self.__params.generate_velocity_noise( self.__params.velocity_perturbations[j], - worker_state.leadtime_since_start_forecast[j] - * self.__config.timestep, + worker_state.leadtime_since_start_forecast[j] * self.__config.timestep, ) ) @@ -2560,9 +2443,7 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( ) # Extrapolate the extrapolation and noise cascade - extrap_kwargs_["displacement_prev"] = ( - worker_state.previous_displacement[j] - ) + extrap_kwargs_["displacement_prev"] = worker_state.previous_displacement[j] extrap_kwargs_noise["displacement_prev"] = ( worker_state.previous_displacement_noise_cascade[j] ) @@ -2608,12 +2489,8 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( worker_state.time_prev_timestep[j] = t + 1 - worker_state.precip_cascades_prev_subtimestep[j] = worker_state.precip_cascades[ - j - ] - worker_state.cascade_noise_prev_subtimestep[j] = ( - worker_state.precip_noise_cascades[j] - ) + worker_state.precip_cascades_prev_subtimestep[j] = worker_state.precip_cascades[j] + worker_state.cascade_noise_prev_subtimestep[j] = worker_state.precip_noise_cascades[j] def __blend_cascades(self, t_sub, j, worker_state): """ @@ -2621,9 +2498,9 @@ def __blend_cascades(self, t_sub, j, worker_state): Computes both full and model-only blends while also blending means and standard deviations across scales. """ - worker_state.subtimestep_index = np.where( - np.array(worker_state.subtimesteps) == t_sub - )[0][0] + worker_state.subtimestep_index = np.where(np.array(worker_state.subtimesteps) == t_sub)[0][ + 0 + ] # First concatenate the cascades and the means and sigmas # precip_models = [n_models,timesteps,n_cascade_levels,m,n] @@ -2634,42 +2511,32 @@ def __blend_cascades(self, t_sub, j, worker_state): ): cascade_stack_all_components = np.concatenate( ( - worker_state.precip_extrapolated_decomp[ - None, worker_state.subtimestep_index - ], + worker_state.precip_extrapolated_decomp[None, worker_state.subtimestep_index], worker_state.precip_models_cascades_timestep, ), axis=0, ) # [(extr_field, n_model_fields), n_cascade_levels, ...] - means_stacked = np.concatenate( ( - worker_state.mean_nowcast_timestep, + worker_state.mean_nowcast_timestep[t_sub, j], worker_state.mean_models_timestep, ), axis=0, ) sigmas_stacked = np.concatenate( ( - worker_state.std_nowcast_timestep, + worker_state.std_nowcast_timestep[t_sub, j], worker_state.std_models_timestep, ), axis=0, ) - elif ( - self.__config.blend_nwp_members - and self.__config.nowcasting_method == "steps" - ): + elif self.__config.blend_nwp_members and self.__config.nowcasting_method == "steps": cascade_stack_all_components = np.concatenate( ( - worker_state.precip_extrapolated_decomp[ - None, worker_state.subtimestep_index - ], + worker_state.precip_extrapolated_decomp[None, worker_state.subtimestep_index], worker_state.precip_models_cascades_timestep, - worker_state.noise_extrapolated_decomp[ - None, worker_state.subtimestep_index - ], + worker_state.noise_extrapolated_decomp[None, worker_state.subtimestep_index], ), axis=0, ) # [(extr_field, n_model_fields, noise), n_cascade_levels, ...] @@ -2691,9 +2558,7 @@ def __blend_cascades(self, t_sub, j, worker_state): elif self.__config.nowcasting_method == "external_nowcast": cascade_stack_all_components = np.concatenate( ( - worker_state.precip_extrapolated_decomp[ - None, worker_state.subtimestep_index - ], + worker_state.precip_extrapolated_decomp[None, worker_state.subtimestep_index], worker_state.precip_models_cascades_timestep[None, j], ), axis=0, @@ -2701,7 +2566,7 @@ def __blend_cascades(self, t_sub, j, worker_state): means_stacked = np.concatenate( ( # worker_state.mean_nowcast_timestep, - worker_state.mean_extrapolation[None, :], + worker_state.mean_nowcast_timestep[t_sub, j], worker_state.mean_models_timestep[None, j], ), axis=0, @@ -2709,7 +2574,7 @@ def __blend_cascades(self, t_sub, j, worker_state): sigmas_stacked = np.concatenate( ( # worker_state.std_nowcast_timestep, - worker_state.std_extrapolation[None, :], + worker_state.std_nowcast_timestep[t_sub, j], worker_state.std_models_timestep[None, j], ), axis=0, @@ -2718,13 +2583,9 @@ def __blend_cascades(self, t_sub, j, worker_state): else: cascade_stack_all_components = np.concatenate( ( - worker_state.precip_extrapolated_decomp[ - None, worker_state.subtimestep_index - ], + worker_state.precip_extrapolated_decomp[None, worker_state.subtimestep_index], worker_state.precip_models_cascades_timestep[None, j], - worker_state.noise_extrapolated_decomp[ - None, worker_state.subtimestep_index - ], + worker_state.noise_extrapolated_decomp[None, worker_state.subtimestep_index], ), axis=0, ) # [(extr_field, n_model_fields, noise), n_cascade_levels, ...] @@ -2773,46 +2634,36 @@ def __blend_cascades(self, t_sub, j, worker_state): # Create weights_with_noise to ensure there is always a 3D weights field, even # if self.__config.nowcasting_method is "external_nowcast" and n_ens_members is 1. worker_state.weights_with_noise = worker_state.weights.copy() - worker_state.weights_model_only_with_noise = ( - worker_state.weights_model_only.copy() - ) + worker_state.weights_model_only_with_noise = worker_state.weights_model_only.copy() if self.__config.nowcasting_method == "external_nowcast": # First determine the weights without noise worker_state.weights = worker_state.weights[:-1, :] / np.sum( worker_state.weights[:-1, :], axis=0 ) - worker_state.weights_model_only = worker_state.weights_model_only[ - :-1, : - ] / np.sum(worker_state.weights_model_only[:-1, :], axis=0) + worker_state.weights_model_only = worker_state.weights_model_only[:-1, :] / np.sum( + worker_state.weights_model_only[:-1, :], axis=0 + ) # Blend the extrapolation, (NWP) model(s) and noise cascades - worker_state.final_blended_forecast_cascades = ( - blending.utils.blend_cascades( - cascades_norm=cascade_stack_all_components, - weights=worker_state.weights, - ) + worker_state.final_blended_forecast_cascades = blending.utils.blend_cascades( + cascades_norm=cascade_stack_all_components, + weights=worker_state.weights, ) # Also blend the cascade without the extrapolation component - worker_state.final_blended_forecast_cascades_mod_only = ( - blending.utils.blend_cascades( - cascades_norm=cascade_stack_all_components[1:, :], - weights=worker_state.weights_model_only, - ) + worker_state.final_blended_forecast_cascades_mod_only = blending.utils.blend_cascades( + cascades_norm=cascade_stack_all_components[1:, :], + weights=worker_state.weights_model_only, ) else: # Blend the extrapolation, (NWP) model(s) and noise cascades - worker_state.final_blended_forecast_cascades = ( - blending.utils.blend_cascades( - cascades_norm=cascade_stack_all_components, - weights=worker_state.weights_with_noise, - ) + worker_state.final_blended_forecast_cascades = blending.utils.blend_cascades( + cascades_norm=cascade_stack_all_components, + weights=worker_state.weights_with_noise, ) # Also blend the cascade without the extrapolation component - worker_state.final_blended_forecast_cascades_mod_only = ( - blending.utils.blend_cascades( - cascades_norm=cascade_stack_all_components[1:, :], - weights=worker_state.weights_model_only, - ) + worker_state.final_blended_forecast_cascades_mod_only = blending.utils.blend_cascades( + cascades_norm=cascade_stack_all_components[1:, :], + weights=worker_state.weights_model_only, ) # Blend the means and standard deviations @@ -2841,36 +2692,28 @@ def __recompose_cascade_to_rainfall_field(self, j, worker_state): means and standard deviations. If using the spectral domain, apply inverse FFT for reconstruction. """ - worker_state.final_blended_forecast_recomposed = ( - blending.utils.recompose_cascade( - combined_cascade=worker_state.final_blended_forecast_cascades, - combined_mean=worker_state.final_blended_forecast_means, - combined_sigma=worker_state.final_blended_forecast_stds, - ) + worker_state.final_blended_forecast_recomposed = blending.utils.recompose_cascade( + combined_cascade=worker_state.final_blended_forecast_cascades, + combined_mean=worker_state.final_blended_forecast_means, + combined_sigma=worker_state.final_blended_forecast_stds, ) # The recomposed cascade without the extrapolation (for NaN filling # outside the radar domain) - worker_state.final_blended_forecast_recomposed_mod_only = ( - blending.utils.recompose_cascade( - combined_cascade=worker_state.final_blended_forecast_cascades_mod_only, - combined_mean=worker_state.final_blended_forecast_means_mod_only, - combined_sigma=worker_state.final_blended_forecast_stds_mod_only, - ) + worker_state.final_blended_forecast_recomposed_mod_only = blending.utils.recompose_cascade( + combined_cascade=worker_state.final_blended_forecast_cascades_mod_only, + combined_mean=worker_state.final_blended_forecast_means_mod_only, + combined_sigma=worker_state.final_blended_forecast_stds_mod_only, ) if self.__config.domain == "spectral": # TODO: Check this! (Only tested with domain == 'spatial') - worker_state.final_blended_forecast_recomposed = self.__params.fft_objs[ - j - ].irfft2(worker_state.final_blended_forecast_recomposed) - worker_state.final_blended_forecast_recomposed_mod_only = ( - self.__params.fft_objs[j].irfft2( - worker_state.final_blended_forecast_recomposed_mod_only - ) + worker_state.final_blended_forecast_recomposed = self.__params.fft_objs[j].irfft2( + worker_state.final_blended_forecast_recomposed ) + worker_state.final_blended_forecast_recomposed_mod_only = self.__params.fft_objs[ + j + ].irfft2(worker_state.final_blended_forecast_recomposed_mod_only) - def __post_process_output( - self, j, t_sub, final_blended_forecast_single_member, worker_state - ): + def __post_process_output(self, j, t_sub, final_blended_forecast_single_member, worker_state): """ Apply post-processing steps to refine the final blended forecast. This involves masking, filling missing data with the blended NWP forecast, @@ -2910,12 +2753,11 @@ def __post_process_output( ) # And the weights for outside the radar domain - weights_probability_matching_mod_only = ( - worker_state.weights_model_only_with_noise[:-1, 1] - ) # Weights without noise, level 2 + weights_probability_matching_mod_only = worker_state.weights_model_only_with_noise[ + :-1, 1 + ] # Weights without noise, level 2 weights_probability_matching_normalized_mod_only = ( - weights_probability_matching_mod_only - / np.sum(weights_probability_matching_mod_only) + weights_probability_matching_mod_only / np.sum(weights_probability_matching_mod_only) ) # Stack the fields if self.__config.blend_nwp_members: @@ -3028,9 +2870,7 @@ def __post_process_output( # apply the precipitation mask to prevent generation of new # precipitation into areas where it was not originally # observed - precip_forecast_min_value = ( - worker_state.final_blended_forecast_recomposed.min() - ) + precip_forecast_min_value = worker_state.final_blended_forecast_recomposed.min() if self.__config.mask_method == "incremental": # The incremental mask is slightly different from the implementation in # nowcasts.steps.py, as it is not computed in the Lagrangian space. Instead, @@ -3038,8 +2878,7 @@ def __post_process_output( # the time step until mask_rim_max. This ensures that for the first t time # steps, the buffer mask keeps increasing. precip_field_mask = ( - precip_forecast_probability_matching_blended - >= self.__params.precip_threshold + precip_forecast_probability_matching_blended >= self.__params.precip_threshold ) # Buffer the mask @@ -3056,9 +2895,7 @@ def __post_process_output( accumulated_mask = dilated_mask.astype(float) # Iteratively dilate the mask and accumulate the results to create a grayscale rim - mask_rim_temp = min( - self.__params.mask_rim + t_sub - 1, self.__params.max_mask_rim - ) + mask_rim_temp = min(self.__params.mask_rim + t_sub - 1, self.__params.max_mask_rim) for _ in range(mask_rim_temp): dilated_mask = binary_dilation(dilated_mask, struct_element) accumulated_mask += dilated_mask @@ -3068,22 +2905,17 @@ def __post_process_output( # Get the final mask worker_state.final_blended_forecast_recomposed = ( precip_forecast_min_value - + ( - worker_state.final_blended_forecast_recomposed - - precip_forecast_min_value - ) + + (worker_state.final_blended_forecast_recomposed - precip_forecast_min_value) * precip_field_mask ) precip_field_mask_temp = ( - worker_state.final_blended_forecast_recomposed - > precip_forecast_min_value + worker_state.final_blended_forecast_recomposed > precip_forecast_min_value ) elif self.__config.mask_method == "obs": # The mask equals the most recent benchmark # rainfall field precip_field_mask_temp = ( - precip_forecast_probability_matching_blended - >= self.__params.precip_threshold + precip_forecast_probability_matching_blended >= self.__params.precip_threshold ) # Set to min value outside of mask @@ -3094,23 +2926,18 @@ def __post_process_output( # If probmatching_method is not None, resample the distribution from # both the extrapolation cascade and the model (NWP) cascade and use # that for the probability matching. - if ( - self.__config.probmatching_method is not None - and self.__config.resample_distribution - ): + if self.__config.probmatching_method is not None and self.__config.resample_distribution: arr1 = worker_state.precip_extrapolated_probability_matching[ worker_state.subtimestep_index ] arr2 = worker_state.precip_models_timestep[j] # resample weights based on cascade level 2. # Areas where one of the fields is nan are not included. - precip_forecast_probability_matching_resampled = ( - probmatching.resample_distributions( - first_array=arr1, - second_array=arr2, - probability_first_array=weights_probability_matching_normalized[0], - randgen=worker_state.randgen_probmatching[j], - ) + precip_forecast_probability_matching_resampled = probmatching.resample_distributions( + first_array=arr1, + second_array=arr2, + probability_first_array=weights_probability_matching_normalized[0], + randgen=worker_state.randgen_probmatching[j], ) else: precip_forecast_probability_matching_resampled = ( @@ -3140,13 +2967,11 @@ def __post_process_output( # Use R_pm_blended as benchmark field and mean_probabiltity_matching_forecast = np.mean( precip_forecast_probability_matching_resampled[ - precip_forecast_probability_matching_resampled - >= self.__params.precip_threshold + precip_forecast_probability_matching_resampled >= self.__params.precip_threshold ] ) no_rain_mask = ( - worker_state.final_blended_forecast_recomposed - >= self.__params.precip_threshold + worker_state.final_blended_forecast_recomposed >= self.__params.precip_threshold ) mean_precip_forecast = np.mean( worker_state.final_blended_forecast_recomposed[no_rain_mask] @@ -3158,9 +2983,7 @@ def __post_process_output( ) precip_forecast_probability_matching_resampled = None - final_blended_forecast_single_member.append( - worker_state.final_blended_forecast_recomposed - ) + final_blended_forecast_single_member.append(worker_state.final_blended_forecast_recomposed) return final_blended_forecast_single_member def __measure_time(self, label, start_time): @@ -3263,10 +3086,8 @@ def forecast( velocity_models: array-like Array of shape (n_models,timestep,2,m,n) containing the x- and y-components of the advection field for the (NWP) model field per forecast lead time. - All values are required to be finite. - - To reduce memory usage, this array - can be given as float32. They will then be converted to float64 before computations + All values are required to be finite. To reduce memory usage, this array can + be given as float32. They will then be converted to float64 before computations to minimize loss in precision. timesteps: int or list of floats Number of time steps to forecast or a list of time steps for which the @@ -3712,9 +3533,7 @@ def calculate_weights_spn(correlations, covariance): cov_matrix_inv = inv(cov_matrix) # The component weights are the dot product between cov_matrix_inv and cor_vec weights = np.dot(cov_matrix_inv, correlations) - weights = np.nan_to_num( - weights, copy=True, nan=10e-5, posinf=10e-5, neginf=10e-5 - ) + weights = np.nan_to_num(weights, copy=True, nan=10e-5, posinf=10e-5, neginf=10e-5) weights_dot_correlations = np.dot(weights, correlations) # If the dot product of the weights with the correlations is # larger than 1.0, we assign a weight of 0.0 to the noise (to make @@ -3790,10 +3609,7 @@ def blend_means_sigmas(means, sigmas, weights): sigmas = np.expand_dims(sigmas, axis=sigmas.ndim) # Weight should have one component more (the noise component) than the # means and sigmas. Check this - if ( - weights.shape[0] - means.shape[0] != 1 - or weights.shape[0] - sigmas.shape[0] != 1 - ): + if weights.shape[0] - means.shape[0] != 1 or weights.shape[0] - sigmas.shape[0] != 1: raise ValueError( "The weights array does not have one (noise) component more than mu and sigma" ) From 3cc1023c9e13922b7008e72a9bafdd4af7c77bd5 Mon Sep 17 00:00:00 2001 From: Ruben Imhoff Date: Wed, 17 Sep 2025 17:30:07 +0200 Subject: [PATCH 06/18] refactor: run black --- pysteps/blending/steps.py | 516 ++++++++++++++++++++++++++------------ 1 file changed, 353 insertions(+), 163 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 26869b12b..2013aef8b 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -572,7 +572,9 @@ def compute_forecast(self): self.__prepare_forecast_loop() self.__initialize_noise_cascades() if self.__config.measure_time: - self.__init_time = self.__measure_time("initialization", self.__start_time_init) + self.__init_time = self.__measure_time( + "initialization", self.__start_time_init + ) self.__blended_nowcast_main_loop() # Stack and return the forecast output @@ -607,10 +609,16 @@ def __blended_nowcast_main_loop(self): starttime_mainloop = time.time() self.__state.extrapolation_kwargs["return_displacement"] = True - self.__state.precip_cascades_prev_subtimestep = deepcopy(self.__state.precip_cascades) - self.__state.cascade_noise_prev_subtimestep = deepcopy(self.__state.precip_noise_cascades) + self.__state.precip_cascades_prev_subtimestep = deepcopy( + self.__state.precip_cascades + ) + self.__state.cascade_noise_prev_subtimestep = deepcopy( + self.__state.precip_noise_cascades + ) - self.__state.time_prev_timestep = [0.0 for j in range(self.__config.n_ens_members)] + self.__state.time_prev_timestep = [ + 0.0 for j in range(self.__config.n_ens_members) + ] self.__state.leadtime_since_start_forecast = [ 0.0 for j in range(self.__config.n_ens_members) ] @@ -643,11 +651,13 @@ def worker(j): if t_sub > 0: self.__blend_cascades(t_sub, j, worker_state) self.__recompose_cascade_to_rainfall_field(j, worker_state) - final_blended_forecast_single_member = self.__post_process_output( - j, - t_sub, - final_blended_forecast_single_member, - worker_state, + final_blended_forecast_single_member = ( + self.__post_process_output( + j, + t_sub, + final_blended_forecast_single_member, + worker_state, + ) ) final_blended_forecast_all_members_one_timestep[j] = ( final_blended_forecast_single_member @@ -675,7 +685,9 @@ def worker(j): print("done.") if self.__config.callback is not None: - precip_forecast_final = np.stack(final_blended_forecast_all_members_one_timestep) + precip_forecast_final = np.stack( + final_blended_forecast_all_members_one_timestep + ) if precip_forecast_final.shape[1] > 0: self.__config.callback(precip_forecast_final.squeeze()) @@ -713,14 +725,19 @@ def __check_inputs(self): raise KeyError( "if precip_nowcast is not None, nowcasting_method should be set to 'external_nowcast' " ) - if self.__config.nowcasting_method == "external_nowcast" and self.__precip_nowcast is None: + if ( + self.__config.nowcasting_method == "external_nowcast" + and self.__precip_nowcast is None + ): raise KeyError( "if nowcasting_method is set to 'external_nowcast', an external precip_nowcast should be provided as variable." ) # Check dimensions of velocity if self.__velocity.ndim != 3: - raise ValueError("velocity must be a three-dimensional array of shape (2, m, n)") + raise ValueError( + "velocity must be a three-dimensional array of shape (2, m, n)" + ) if self.__velocity_models.ndim != 5: raise ValueError( "velocity_models must be a five-dimensional array of shape (n_models, timestep, 2, m, n)" @@ -745,13 +762,17 @@ def __check_inputs(self): if isinstance(self.__timesteps, list): self.__params.time_steps_is_list = True if not sorted(self.__timesteps) == self.__timesteps: - raise ValueError("timesteps is not in ascending order", self.__timesteps) + raise ValueError( + "timesteps is not in ascending order", self.__timesteps + ) if self.__precip_models.shape[1] != math.ceil(self.__timesteps[-1]) + 1: raise ValueError( "precip_models does not contain sufficient lead times for this forecast" ) self.__params.original_timesteps = [0] + list(self.__timesteps) - self.__timesteps = nowcast_utils.binned_timesteps(self.__params.original_timesteps) + self.__timesteps = nowcast_utils.binned_timesteps( + self.__params.original_timesteps + ) else: self.__params.time_steps_is_list = False if self.__precip_models.shape[1] != self.__timesteps + 1: @@ -784,7 +805,9 @@ def __check_inputs(self): if self.__config.extrapolation_kwargs is None: self.__state.extrapolation_kwargs = dict() else: - self.__state.extrapolation_kwargs = deepcopy(self.__config.extrapolation_kwargs) + self.__state.extrapolation_kwargs = deepcopy( + self.__config.extrapolation_kwargs + ) if self.__config.filter_kwargs is None: self.__params.filter_kwargs = dict() @@ -805,9 +828,13 @@ def __check_inputs(self): if self.__config.climatology_kwargs is None: # Make sure clim_kwargs at least contains the number of models - self.__params.climatology_kwargs = dict({"n_models": self.__precip_models.shape[0]}) + self.__params.climatology_kwargs = dict( + {"n_models": self.__precip_models.shape[0]} + ) else: - self.__params.climatology_kwargs = deepcopy(self.__config.climatology_kwargs) + self.__params.climatology_kwargs = deepcopy( + self.__config.climatology_kwargs + ) if self.__config.mask_kwargs is None: self.__params.mask_kwargs = dict() @@ -828,7 +855,10 @@ def __check_inputs(self): if self.__config.conditional and self.__params.precip_threshold is None: raise ValueError("conditional=True but precip_thr is not set") - if self.__config.mask_method is not None and self.__params.precip_threshold is None: + if ( + self.__config.mask_method is not None + and self.__params.precip_threshold is None + ): raise ValueError("mask_method!=None but precip_thr=None") if self.__config.noise_stddev_adj not in ["auto", "fixed", None]: @@ -839,13 +869,17 @@ def __check_inputs(self): if self.__config.kmperpixel is None: if self.__config.velocity_perturbation_method is not None: - raise ValueError("velocity_perturbation_method is set but kmperpixel=None") + raise ValueError( + "velocity_perturbation_method is set but kmperpixel=None" + ) if self.__config.mask_method == "incremental": raise ValueError("mask_method='incremental' but kmperpixel=None") if self.__config.timestep is None: if self.__config.velocity_perturbation_method is not None: - raise ValueError("velocity_perturbation_method is set but timestep=None") + raise ValueError( + "velocity_perturbation_method is set but timestep=None" + ) if self.__config.mask_method == "incremental": raise ValueError("mask_method='incremental' but timestep=None") @@ -860,7 +894,9 @@ def __print_forecast_info(self): print("Inputs") print("------") print(f"forecast issue time: {self.__issuetime.isoformat()}") - print(f"input dimensions: {self.__precip.shape[1]}x{self.__precip.shape[2]}") + print( + f"input dimensions: {self.__precip.shape[1]}x{self.__precip.shape[2]}" + ) if self.__precip_nowcast is not None: print( f"input dimensions pre-computed nowcast: {self.__precip_nowcast.shape[2]}x{self.__precip_nowcast.shape[3]}" @@ -875,7 +911,9 @@ def __print_forecast_info(self): print("-----------------------") print(f"number of (NWP) models: {self.__precip_models.shape[0]}") print(f"blend (NWP) model members: {self.__config.blend_nwp_members}") - print(f"decompose (NWP) models: {'yes' if self.__precip_models.ndim == 4 else 'no'}") + print( + f"decompose (NWP) models: {'yes' if self.__precip_models.ndim == 4 else 'no'}" + ) print("") print("Methods") @@ -885,10 +923,16 @@ def __print_forecast_info(self): print(f"decomposition: {self.__config.decomposition_method}") print(f"nowcasting algorithm: {self.__config.nowcasting_method}") print(f"noise generator: {self.__config.noise_method}") - print(f"noise adjustment: {'yes' if self.__config.noise_stddev_adj else 'no'}") - print(f"velocity perturbator: {self.__config.velocity_perturbation_method}") + print( + f"noise adjustment: {'yes' if self.__config.noise_stddev_adj else 'no'}" + ) + print( + f"velocity perturbator: {self.__config.velocity_perturbation_method}" + ) print(f"blending weights method: {self.__config.weights_method}") - print(f"conditional statistics: {'yes' if self.__config.conditional else 'no'}") + print( + f"conditional statistics: {'yes' if self.__config.conditional else 'no'}" + ) print(f"precip. mask method: {self.__config.mask_method}") print(f"probability matching: {self.__config.probmatching_method}") print(f"FFT method: {self.__config.fft_method}") @@ -1028,24 +1072,29 @@ def transform_to_lagrangian(precip, i): for i in range(self.__config.ar_order): res.append(dask.delayed(transform_to_lagrangian)(self.__precip, i)) num_workers_ = ( - len(res) if self.__config.num_workers > len(res) else self.__config.num_workers + len(res) + if self.__config.num_workers > len(res) + else self.__config.num_workers ) self.__precip = np.stack( - list(dask.compute(*res, num_workers=num_workers_)) + [self.__precip[-1, :, :]] + list(dask.compute(*res, num_workers=num_workers_)) + + [self.__precip[-1, :, :]] ) # Replace non-finite values with the minimum value for each field self.__precip = self.__precip.copy() for i in range(self.__precip.shape[0]): - self.__precip[i, ~np.isfinite(self.__precip[i, :])] = np.nanmin(self.__precip[i, :]) + self.__precip[i, ~np.isfinite(self.__precip[i, :])] = np.nanmin( + self.__precip[i, :] + ) if self.__precip_nowcast is not None: self.__precip_nowcast = self.__precip_nowcast.copy() # TODO: check this function for ensemble nowcasts -> now made it always work with 4 dimentions for m in range(self.__precip_nowcast.shape[0]): for t in range(self.__precip_nowcast.shape[1]): - self.__precip_nowcast[m, t, ~np.isfinite(self.__precip_nowcast[m, t, :, :])] = ( - np.nanmin(self.__precip_nowcast[m, t, :, :]) - ) + self.__precip_nowcast[ + m, t, ~np.isfinite(self.__precip_nowcast[m, t, :, :]) + ] = np.nanmin(self.__precip_nowcast[m, t, :, :]) # Perform the cascade decomposition for the input precip fields and, # if necessary, for the (NWP) model fields @@ -1082,16 +1131,22 @@ def transform_to_lagrangian(precip, i): ) else: with ThreadPool(self.__config.num_workers) as pool: - precip_nowcast_decomp, precip_nowcast_means, precip_nowcast_stds = pool.map( - partial(self.__decompose_member), - list( - self.__precip_nowcast - ), # each entry = one ensemble member (shape [n_timesteps, x, y]) + precip_nowcast_decomp, precip_nowcast_means, precip_nowcast_stds = ( + pool.map( + partial(self.__decompose_member), + list(self.__precip_nowcast), + ) ) - self.__state.precip_nowcast_cascades = np.array([precip_nowcast_decomp]).swapaxes(1, 2) - self.__state.mean_nowcast_timestep = np.array([precip_nowcast_means]).swapaxes(1, 2) - self.__state.std_nowcast_timestep = np.array([precip_nowcast_stds]).swapaxes(1, 2) + self.__state.precip_nowcast_cascades = np.array( + [precip_nowcast_decomp] + ).swapaxes(1, 2) + self.__state.mean_nowcast_timestep = np.array( + [precip_nowcast_means] + ).swapaxes(1, 2) + self.__state.std_nowcast_timestep = np.array( + [precip_nowcast_stds] + ).swapaxes(1, 2) # If necessary, recompose (NWP) model forecasts self.__state.precip_models_cascades = None @@ -1160,7 +1215,9 @@ def __zero_precipitation_forecast(self): when no precipitation above the threshold is detected. The forecast is optionally returned or passed to a callback. """ - print("No precipitation above the threshold found in both the radar and NWP fields") + print( + "No precipitation above the threshold found in both the radar and NWP fields" + ) print("The resulting forecast will contain only zeros") # Create the output list precip_forecast = [[] for j in range(self.__config.n_ens_members)] @@ -1191,7 +1248,10 @@ def __zero_precipitation_forecast(self): if self.__config.return_output: precip_forecast_all_members_all_times = np.stack( - [np.stack(precip_forecast[j]) for j in range(self.__config.n_ens_members)] + [ + np.stack(precip_forecast[j]) + for j in range(self.__config.n_ens_members) + ] ) if self.__config.measure_time: @@ -1215,7 +1275,9 @@ def __prepare_nowcast_for_zero_radar(self): # forecast. I.e., velocity (radar) equals velocity_model at the first time # step. # Use the velocity from velocity_models at time step 0 - self.__velocity = self.__velocity_models[:, 0, :, :, :].astype(np.float64, copy=False) + self.__velocity = self.__velocity_models[:, 0, :, :, :].astype( + np.float64, copy=False + ) # Take the average over the first axis, which corresponds to n_models # (hence, the model average) self.__velocity = np.mean(self.__velocity, axis=0) @@ -1241,7 +1303,9 @@ def __prepare_nowcast_for_zero_radar(self): max_rain_pixels = rain_pixels max_rain_pixels_j = j max_rain_pixels_t = t - self.__state.precip_noise_input = self.__precip_models[max_rain_pixels_j][max_rain_pixels_t] + self.__state.precip_noise_input = self.__precip_models[max_rain_pixels_j][ + max_rain_pixels_t + ] self.__state.precip_noise_input = self.__state.precip_noise_input.astype( np.float64, copy=False ) @@ -1249,10 +1313,12 @@ def __prepare_nowcast_for_zero_radar(self): # If zero_precip_radar, make sure that precip_cascade does not contain # only nans or infs. If so, fill it with the zero value. if self.__state.precip_models_cascades is not None: - self.__state.precip_cascades[~np.isfinite(self.__state.precip_cascades)] = np.nanmin( - self.__state.precip_models_cascades[max_rain_pixels_j, max_rain_pixels_t][ - "cascade_levels" - ] + self.__state.precip_cascades[~np.isfinite(self.__state.precip_cascades)] = ( + np.nanmin( + self.__state.precip_models_cascades[ + max_rain_pixels_j, max_rain_pixels_t + ]["cascade_levels"] + ) ) else: precip_models_cascade_timestep = self.__params.decomposition_method( @@ -1264,13 +1330,15 @@ def __prepare_nowcast_for_zero_radar(self): compute_stats=True, compact_output=True, )["cascade_levels"] - self.__state.precip_cascades[~np.isfinite(self.__state.precip_cascades)] = np.nanmin( - precip_models_cascade_timestep + self.__state.precip_cascades[~np.isfinite(self.__state.precip_cascades)] = ( + np.nanmin(precip_models_cascade_timestep) ) # Make sure precip_noise_input is three-dimensional if len(self.__state.precip_noise_input.shape) != 3: - self.__state.precip_noise_input = self.__state.precip_noise_input[np.newaxis, :, :] + self.__state.precip_noise_input = self.__state.precip_noise_input[ + np.newaxis, :, : + ] def __initialize_noise(self): """ @@ -1279,7 +1347,9 @@ def __initialize_noise(self): """ if self.__config.noise_method is not None: # get methods for perturbations - init_noise, self.__params.noise_generator = noise.get_method(self.__config.noise_method) + init_noise, self.__params.noise_generator = noise.get_method( + self.__config.noise_method + ) # initialize the perturbation generator for the precipitation field self.__params.perturbation_generator = init_noise( @@ -1391,11 +1461,15 @@ def __estimate_ar_parameters_radar(self): # adjust the lag-2 correlation coefficient to ensure that the AR(p) # process is stationary for i in range(self.__config.n_cascade_levels): - GAMMA[i, 1] = autoregression.adjust_lag2_corrcoef2(GAMMA[i, 0], GAMMA[i, 1]) + GAMMA[i, 1] = autoregression.adjust_lag2_corrcoef2( + GAMMA[i, 0], GAMMA[i, 1] + ) # estimate the parameters of the AR(p) model from the auto-correlation # coefficients - self.__params.PHI = np.empty((self.__config.n_cascade_levels, self.__config.ar_order + 1)) + self.__params.PHI = np.empty( + (self.__config.n_cascade_levels, self.__config.ar_order + 1) + ) for i in range(self.__config.n_cascade_levels): self.__params.PHI[i, :] = autoregression.estimate_ar_params_yw(GAMMA[i, :]) @@ -1486,12 +1560,16 @@ def __prepare_forecast_loop(self): self.__state.previous_displacement_prob_matching = np.stack( [None for j in range(self.__config.n_ens_members)] ) - self.__state.final_blended_forecast = [[] for j in range(self.__config.n_ens_members)] + self.__state.final_blended_forecast = [ + [] for j in range(self.__config.n_ens_members) + ] if self.__config.mask_method == "incremental": # get mask parameters self.__params.mask_rim = self.__params.mask_kwargs.get("mask_rim", 10) - self.__params.max_mask_rim = self.__params.mask_kwargs.get("max_mask_rim", 10) + self.__params.max_mask_rim = self.__params.mask_kwargs.get( + "max_mask_rim", 10 + ) mask_f = self.__params.mask_kwargs.get("mask_f", 1.0) # initialize the structuring element struct = generate_binary_structure(2, 1) @@ -1520,8 +1598,12 @@ def __prepare_forecast_loop(self): # initizalize the current and previous extrapolation forecast scale for the nowcasting component # phi1 / (1 - phi2), see BPS2004 - self.__state.rho_extrap_cascade_prev = np.repeat(1.0, self.__params.PHI.shape[0]) - self.__state.rho_extrap_cascade = self.__params.PHI[:, 0] / (1.0 - self.__params.PHI[:, 1]) + self.__state.rho_extrap_cascade_prev = np.repeat( + 1.0, self.__params.PHI.shape[0] + ) + self.__state.rho_extrap_cascade = self.__params.PHI[:, 0] / ( + 1.0 - self.__params.PHI[:, 1] + ) def __initialize_noise_cascades(self): """ @@ -1529,7 +1611,9 @@ def __initialize_noise_cascades(self): We also need to return the mean and standard deviations of the noise for the recombination of the noise before advecting it. """ - self.__state.precip_noise_cascades = np.zeros(self.__state.precip_cascades.shape) + self.__state.precip_noise_cascades = np.zeros( + self.__state.precip_cascades.shape + ) self.__state.precip_mean_noise = np.zeros( (self.__config.n_ens_members, self.__config.n_cascade_levels) ) @@ -1650,13 +1734,17 @@ def __decompose_nwp_if_needed_and_fill_nans_in_nwp(self, t): self.__state.precip_models_cascades_timestep[ ~np.isfinite(self.__state.precip_models_cascades_timestep) ] = min_cascade - self.__state.precip_models_timestep[~np.isfinite(self.__state.precip_models_timestep)] = ( - min_precip - ) + self.__state.precip_models_timestep[ + ~np.isfinite(self.__state.precip_models_timestep) + ] = min_precip # Also set any nans or infs in the mean and sigma of the cascade to # respectively 0.0 and 1.0 - self.__state.mean_models_timestep[~np.isfinite(self.__state.mean_models_timestep)] = 0.0 - self.__state.std_models_timestep[~np.isfinite(self.__state.std_models_timestep)] = 0.0 + self.__state.mean_models_timestep[ + ~np.isfinite(self.__state.mean_models_timestep) + ] = 0.0 + self.__state.std_models_timestep[ + ~np.isfinite(self.__state.std_models_timestep) + ] = 0.0 def __find_nowcast_NWP_combination(self, t): """ @@ -1668,9 +1756,9 @@ def __find_nowcast_NWP_combination(self, t): np.float64, copy=False ) - self.__state.velocity_models_timestep = self.__velocity_models[:, t, :, :, :].astype( - np.float64, copy=False - ) + self.__state.velocity_models_timestep = self.__velocity_models[ + :, t, :, :, : + ].astype(np.float64, copy=False) # Make sure the number of model members is not larger than or equal to n_ens_members n_model_members = self.__state.precip_models_cascades_timestep.shape[0] if n_model_members > self.__config.n_ens_members: @@ -1698,7 +1786,9 @@ def __find_nowcast_NWP_combination(self, t): # for indexing the right climatological skill file when pysteps calculates # the blended forecast in parallel. if n_model_members > 1: - self.__state.mapping_list_NWP_member_to_ensemble_member = np.arange(n_model_members) + self.__state.mapping_list_NWP_member_to_ensemble_member = np.arange( + n_model_members + ) else: self.__state.mapping_list_NWP_member_to_ensemble_member = [0] @@ -1802,7 +1892,9 @@ def __find_nowcast_NWP_combination(self, t): self.__state.mean_nowcast_timestep = np.repeat( self.__state.mean_nowcast_timestep, repeats, axis=0 ) - self.__state.std_nowcast = np.repeat(self.__state.std_nowcast, repeats, axis=0) + self.__state.std_nowcast = np.repeat( + self.__state.std_nowcast, repeats, axis=0 + ) self.__state.velocity_nowcast = np.repeat( self.__state.velocity_nowcast, repeats, axis=0 ) @@ -1820,7 +1912,9 @@ def __find_nowcast_NWP_combination(self, t): # for indexing the right climatological skill file when pysteps calculates # the blended forecast in parallel. if n_model_members > 1: - self.__state.mapping_list_NWP_member_to_ensemble_member = np.arange(n_model_members) + self.__state.mapping_list_NWP_member_to_ensemble_member = np.arange( + n_model_members + ) else: self.__state.mapping_list_NWP_member_to_ensemble_member = [0] @@ -1898,26 +1992,34 @@ def __determine_skill_for_current_timestep(self, t): if t == 0: # Calculate the initial skill of the (NWP) model forecasts at t=0. self.__params.rho_nwp_models = [] - for model_index in range(self.__state.precip_models_cascades_timestep.shape[0]): + for model_index in range( + self.__state.precip_models_cascades_timestep.shape[0] + ): rho_value = blending.skill_scores.spatial_correlation( obs=self.__state.precip_cascades[0, :, -1, :, :].copy(), - mod=self.__state.precip_models_cascades_timestep[model_index, :, :, :].copy(), + mod=self.__state.precip_models_cascades_timestep[ + model_index, :, :, : + ].copy(), domain_mask=self.__params.domain_mask, ) self.__params.rho_nwp_models.append(rho_value) self.__params.rho_nwp_models = np.stack(self.__params.rho_nwp_models) # Ensure that the model skill decreases with increasing scale level. - for model_index in range(self.__state.precip_models_cascades_timestep.shape[0]): - for i in range(1, self.__state.precip_models_cascades_timestep.shape[1]): + for model_index in range( + self.__state.precip_models_cascades_timestep.shape[0] + ): + for i in range( + 1, self.__state.precip_models_cascades_timestep.shape[1] + ): if ( self.__params.rho_nwp_models[model_index, i] > self.__params.rho_nwp_models[model_index, i - 1] ): # Set it equal to the previous scale level - self.__params.rho_nwp_models[model_index, i] = self.__params.rho_nwp_models[ - model_index, i - 1 - ] + self.__params.rho_nwp_models[model_index, i] = ( + self.__params.rho_nwp_models[model_index, i - 1] + ) # Save this in the climatological skill file blending.clim.save_skill( @@ -1986,7 +2088,9 @@ def __determine_weights_per_component(self, worker_state): # selected, weights will be overwritten with those weights prior to # blending step. # weight = [(extr_field, n_model_fields, noise), n_cascade_levels, ...] - worker_state.weights = calculate_weights_bps(worker_state.rho_final_blended_forecast) + worker_state.weights = calculate_weights_bps( + worker_state.rho_final_blended_forecast + ) # The model only weights if self.__config.weights_method == "bps": @@ -2019,7 +2123,9 @@ def __determine_weights_per_component(self, worker_state): n_model, i, :, : ].flatten() for n_model in range( - worker_state.precip_models_cascades_timestep.shape[0] + worker_state.precip_models_cascades_timestep.shape[ + 0 + ] ) ] ) @@ -2036,7 +2142,8 @@ def __determine_weights_per_component(self, worker_state): ) else: raise ValueError( - "Unknown weights method %s: must be 'bps' or 'spn'" % self.__config.weights_method + "Unknown weights method %s: must be 'bps' or 'spn'" + % self.__config.weights_method ) def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): @@ -2078,7 +2185,9 @@ def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): print("Using nowcasting method:", self.__config.nowcasting_method) for i in range(self.__config.n_cascade_levels): # Use a deterministic Externally computed nowcasting model - worker_state.precip_cascades[j][i] = self.__state.precip_nowcast_cascades[j][i][t] + worker_state.precip_cascades[j][i] = ( + self.__state.precip_nowcast_cascades[j][i][t] + ) # Follow the 'standard' STEPS blending approach as described in :cite:`Imhoff2023` elif self.__config.nowcasting_method == "steps": @@ -2089,8 +2198,10 @@ def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): epsilon_decomposed is not None or self.__config.velocity_perturbation_method is not None ): - worker_state.precip_cascades[j][i] = autoregression.iterate_ar_model( - worker_state.precip_cascades[j][i], self.__params.PHI[i, :] + worker_state.precip_cascades[j][i] = ( + autoregression.iterate_ar_model( + worker_state.precip_cascades[j][i], self.__params.PHI[i, :] + ) ) # Renormalize the cascade worker_state.precip_cascades[j][i][1] /= np.std( @@ -2115,10 +2226,12 @@ def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): epsilon_temp = None # apply AR(p) process to noise cascade level # (Returns zero noise if epsilon_decomposed is None) - worker_state.precip_noise_cascades[j][i] = autoregression.iterate_ar_model( - worker_state.precip_noise_cascades[j][i], - self.__params.PHI[i, :], - eps=epsilon_temp, + worker_state.precip_noise_cascades[j][i] = ( + autoregression.iterate_ar_model( + worker_state.precip_noise_cascades[j][i], + self.__params.PHI[i, :], + eps=epsilon_temp, + ) ) epsilon_decomposed = None @@ -2147,7 +2260,9 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( precip_extrapolated_decomp = worker_state.precip_cascades[j][i][-1, :] worker_state.time_prev_timestep[j] = t + 1 - worker_state.precip_extrapolated_decomp.append(precip_extrapolated_decomp.copy()) + worker_state.precip_extrapolated_decomp.append( + precip_extrapolated_decomp.copy() + ) worker_state.precip_extrapolated_probability_matching.append( precip_extrapolated_decomp.copy() ) @@ -2197,7 +2312,9 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( noise_cascade_subtimestep = np.stack(noise_cascade_subtimestep) t_diff_prev_subtimestep = t_sub - worker_state.time_prev_timestep[j] - worker_state.leadtime_since_start_forecast[j] += t_diff_prev_subtimestep + worker_state.leadtime_since_start_forecast[ + j + ] += t_diff_prev_subtimestep # compute the perturbed motion field - include the NWP # velocities and the weights. Note that we only perturb @@ -2247,10 +2364,12 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( # This is needed to remove the interpolation artifacts. # In addition, the number of extrapolations is greatly reduced # A. Radar Rain - precip_forecast_recomp_subtimestep = blending.utils.recompose_cascade( - combined_cascade=precip_forecast_cascade_subtimestep, - combined_mean=worker_state.mean_extrapolation, - combined_sigma=worker_state.std_extrapolation, + precip_forecast_recomp_subtimestep = ( + blending.utils.recompose_cascade( + combined_cascade=precip_forecast_cascade_subtimestep, + combined_mean=worker_state.mean_extrapolation, + combined_sigma=worker_state.std_extrapolation, + ) ) # Make sure we have values outside the mask if self.__params.zero_precip_radar: @@ -2262,7 +2381,9 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( neginf=self.__params.precip_zerovalue, ) # Put back the mask - precip_forecast_recomp_subtimestep[self.__params.domain_mask] = np.nan + precip_forecast_recomp_subtimestep[self.__params.domain_mask] = ( + np.nan + ) worker_state.extrapolation_kwargs["displacement_prev"] = ( worker_state.previous_displacement[j] ) @@ -2337,13 +2458,17 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( compact_output=True, )["cascade_levels"] for i in range(self.__config.n_cascade_levels): - noise_extrapolated_decomp[i] *= self.__params.noise_std_coeffs[i] + noise_extrapolated_decomp[i] *= self.__params.noise_std_coeffs[ + i + ] # Append the results to the output lists worker_state.precip_extrapolated_decomp.append( precip_extrapolated_decomp.copy() ) - worker_state.noise_extrapolated_decomp.append(noise_extrapolated_decomp.copy()) + worker_state.noise_extrapolated_decomp.append( + noise_extrapolated_decomp.copy() + ) precip_forecast_cascade_subtimestep = None precip_forecast_recomp_subtimestep = None precip_forecast_extrapolated_recomp_subtimestep_temp = None @@ -2364,9 +2489,9 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( ) # Apply the domain mask to the extrapolation component precip_forecast_temp_for_probability_matching = self.__precip.copy() - precip_forecast_temp_for_probability_matching[self.__params.domain_mask] = ( - np.nan - ) + precip_forecast_temp_for_probability_matching[ + self.__params.domain_mask + ] = np.nan ( precip_forecast_extrapolated_probability_matching_temp, @@ -2409,7 +2534,8 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( self.__velocity + self.__params.generate_velocity_noise( self.__params.velocity_perturbations[j], - worker_state.leadtime_since_start_forecast[j] * self.__config.timestep, + worker_state.leadtime_since_start_forecast[j] + * self.__config.timestep, ) ) @@ -2443,7 +2569,9 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( ) # Extrapolate the extrapolation and noise cascade - extrap_kwargs_["displacement_prev"] = worker_state.previous_displacement[j] + extrap_kwargs_["displacement_prev"] = ( + worker_state.previous_displacement[j] + ) extrap_kwargs_noise["displacement_prev"] = ( worker_state.previous_displacement_noise_cascade[j] ) @@ -2489,8 +2617,12 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( worker_state.time_prev_timestep[j] = t + 1 - worker_state.precip_cascades_prev_subtimestep[j] = worker_state.precip_cascades[j] - worker_state.cascade_noise_prev_subtimestep[j] = worker_state.precip_noise_cascades[j] + worker_state.precip_cascades_prev_subtimestep[j] = worker_state.precip_cascades[ + j + ] + worker_state.cascade_noise_prev_subtimestep[j] = ( + worker_state.precip_noise_cascades[j] + ) def __blend_cascades(self, t_sub, j, worker_state): """ @@ -2498,9 +2630,9 @@ def __blend_cascades(self, t_sub, j, worker_state): Computes both full and model-only blends while also blending means and standard deviations across scales. """ - worker_state.subtimestep_index = np.where(np.array(worker_state.subtimesteps) == t_sub)[0][ - 0 - ] + worker_state.subtimestep_index = np.where( + np.array(worker_state.subtimesteps) == t_sub + )[0][0] # First concatenate the cascades and the means and sigmas # precip_models = [n_models,timesteps,n_cascade_levels,m,n] @@ -2511,7 +2643,9 @@ def __blend_cascades(self, t_sub, j, worker_state): ): cascade_stack_all_components = np.concatenate( ( - worker_state.precip_extrapolated_decomp[None, worker_state.subtimestep_index], + worker_state.precip_extrapolated_decomp[ + None, worker_state.subtimestep_index + ], worker_state.precip_models_cascades_timestep, ), axis=0, @@ -2531,12 +2665,19 @@ def __blend_cascades(self, t_sub, j, worker_state): axis=0, ) - elif self.__config.blend_nwp_members and self.__config.nowcasting_method == "steps": + elif ( + self.__config.blend_nwp_members + and self.__config.nowcasting_method == "steps" + ): cascade_stack_all_components = np.concatenate( ( - worker_state.precip_extrapolated_decomp[None, worker_state.subtimestep_index], + worker_state.precip_extrapolated_decomp[ + None, worker_state.subtimestep_index + ], worker_state.precip_models_cascades_timestep, - worker_state.noise_extrapolated_decomp[None, worker_state.subtimestep_index], + worker_state.noise_extrapolated_decomp[ + None, worker_state.subtimestep_index + ], ), axis=0, ) # [(extr_field, n_model_fields, noise), n_cascade_levels, ...] @@ -2558,7 +2699,9 @@ def __blend_cascades(self, t_sub, j, worker_state): elif self.__config.nowcasting_method == "external_nowcast": cascade_stack_all_components = np.concatenate( ( - worker_state.precip_extrapolated_decomp[None, worker_state.subtimestep_index], + worker_state.precip_extrapolated_decomp[ + None, worker_state.subtimestep_index + ], worker_state.precip_models_cascades_timestep[None, j], ), axis=0, @@ -2583,9 +2726,13 @@ def __blend_cascades(self, t_sub, j, worker_state): else: cascade_stack_all_components = np.concatenate( ( - worker_state.precip_extrapolated_decomp[None, worker_state.subtimestep_index], + worker_state.precip_extrapolated_decomp[ + None, worker_state.subtimestep_index + ], worker_state.precip_models_cascades_timestep[None, j], - worker_state.noise_extrapolated_decomp[None, worker_state.subtimestep_index], + worker_state.noise_extrapolated_decomp[ + None, worker_state.subtimestep_index + ], ), axis=0, ) # [(extr_field, n_model_fields, noise), n_cascade_levels, ...] @@ -2634,36 +2781,46 @@ def __blend_cascades(self, t_sub, j, worker_state): # Create weights_with_noise to ensure there is always a 3D weights field, even # if self.__config.nowcasting_method is "external_nowcast" and n_ens_members is 1. worker_state.weights_with_noise = worker_state.weights.copy() - worker_state.weights_model_only_with_noise = worker_state.weights_model_only.copy() + worker_state.weights_model_only_with_noise = ( + worker_state.weights_model_only.copy() + ) if self.__config.nowcasting_method == "external_nowcast": # First determine the weights without noise worker_state.weights = worker_state.weights[:-1, :] / np.sum( worker_state.weights[:-1, :], axis=0 ) - worker_state.weights_model_only = worker_state.weights_model_only[:-1, :] / np.sum( - worker_state.weights_model_only[:-1, :], axis=0 - ) + worker_state.weights_model_only = worker_state.weights_model_only[ + :-1, : + ] / np.sum(worker_state.weights_model_only[:-1, :], axis=0) # Blend the extrapolation, (NWP) model(s) and noise cascades - worker_state.final_blended_forecast_cascades = blending.utils.blend_cascades( - cascades_norm=cascade_stack_all_components, - weights=worker_state.weights, + worker_state.final_blended_forecast_cascades = ( + blending.utils.blend_cascades( + cascades_norm=cascade_stack_all_components, + weights=worker_state.weights, + ) ) # Also blend the cascade without the extrapolation component - worker_state.final_blended_forecast_cascades_mod_only = blending.utils.blend_cascades( - cascades_norm=cascade_stack_all_components[1:, :], - weights=worker_state.weights_model_only, + worker_state.final_blended_forecast_cascades_mod_only = ( + blending.utils.blend_cascades( + cascades_norm=cascade_stack_all_components[1:, :], + weights=worker_state.weights_model_only, + ) ) else: # Blend the extrapolation, (NWP) model(s) and noise cascades - worker_state.final_blended_forecast_cascades = blending.utils.blend_cascades( - cascades_norm=cascade_stack_all_components, - weights=worker_state.weights_with_noise, + worker_state.final_blended_forecast_cascades = ( + blending.utils.blend_cascades( + cascades_norm=cascade_stack_all_components, + weights=worker_state.weights_with_noise, + ) ) # Also blend the cascade without the extrapolation component - worker_state.final_blended_forecast_cascades_mod_only = blending.utils.blend_cascades( - cascades_norm=cascade_stack_all_components[1:, :], - weights=worker_state.weights_model_only, + worker_state.final_blended_forecast_cascades_mod_only = ( + blending.utils.blend_cascades( + cascades_norm=cascade_stack_all_components[1:, :], + weights=worker_state.weights_model_only, + ) ) # Blend the means and standard deviations @@ -2692,28 +2849,36 @@ def __recompose_cascade_to_rainfall_field(self, j, worker_state): means and standard deviations. If using the spectral domain, apply inverse FFT for reconstruction. """ - worker_state.final_blended_forecast_recomposed = blending.utils.recompose_cascade( - combined_cascade=worker_state.final_blended_forecast_cascades, - combined_mean=worker_state.final_blended_forecast_means, - combined_sigma=worker_state.final_blended_forecast_stds, + worker_state.final_blended_forecast_recomposed = ( + blending.utils.recompose_cascade( + combined_cascade=worker_state.final_blended_forecast_cascades, + combined_mean=worker_state.final_blended_forecast_means, + combined_sigma=worker_state.final_blended_forecast_stds, + ) ) # The recomposed cascade without the extrapolation (for NaN filling # outside the radar domain) - worker_state.final_blended_forecast_recomposed_mod_only = blending.utils.recompose_cascade( - combined_cascade=worker_state.final_blended_forecast_cascades_mod_only, - combined_mean=worker_state.final_blended_forecast_means_mod_only, - combined_sigma=worker_state.final_blended_forecast_stds_mod_only, + worker_state.final_blended_forecast_recomposed_mod_only = ( + blending.utils.recompose_cascade( + combined_cascade=worker_state.final_blended_forecast_cascades_mod_only, + combined_mean=worker_state.final_blended_forecast_means_mod_only, + combined_sigma=worker_state.final_blended_forecast_stds_mod_only, + ) ) if self.__config.domain == "spectral": # TODO: Check this! (Only tested with domain == 'spatial') - worker_state.final_blended_forecast_recomposed = self.__params.fft_objs[j].irfft2( - worker_state.final_blended_forecast_recomposed - ) - worker_state.final_blended_forecast_recomposed_mod_only = self.__params.fft_objs[ + worker_state.final_blended_forecast_recomposed = self.__params.fft_objs[ j - ].irfft2(worker_state.final_blended_forecast_recomposed_mod_only) + ].irfft2(worker_state.final_blended_forecast_recomposed) + worker_state.final_blended_forecast_recomposed_mod_only = ( + self.__params.fft_objs[j].irfft2( + worker_state.final_blended_forecast_recomposed_mod_only + ) + ) - def __post_process_output(self, j, t_sub, final_blended_forecast_single_member, worker_state): + def __post_process_output( + self, j, t_sub, final_blended_forecast_single_member, worker_state + ): """ Apply post-processing steps to refine the final blended forecast. This involves masking, filling missing data with the blended NWP forecast, @@ -2753,11 +2918,12 @@ def __post_process_output(self, j, t_sub, final_blended_forecast_single_member, ) # And the weights for outside the radar domain - weights_probability_matching_mod_only = worker_state.weights_model_only_with_noise[ - :-1, 1 - ] # Weights without noise, level 2 + weights_probability_matching_mod_only = ( + worker_state.weights_model_only_with_noise[:-1, 1] + ) # Weights without noise, level 2 weights_probability_matching_normalized_mod_only = ( - weights_probability_matching_mod_only / np.sum(weights_probability_matching_mod_only) + weights_probability_matching_mod_only + / np.sum(weights_probability_matching_mod_only) ) # Stack the fields if self.__config.blend_nwp_members: @@ -2870,7 +3036,9 @@ def __post_process_output(self, j, t_sub, final_blended_forecast_single_member, # apply the precipitation mask to prevent generation of new # precipitation into areas where it was not originally # observed - precip_forecast_min_value = worker_state.final_blended_forecast_recomposed.min() + precip_forecast_min_value = ( + worker_state.final_blended_forecast_recomposed.min() + ) if self.__config.mask_method == "incremental": # The incremental mask is slightly different from the implementation in # nowcasts.steps.py, as it is not computed in the Lagrangian space. Instead, @@ -2878,7 +3046,8 @@ def __post_process_output(self, j, t_sub, final_blended_forecast_single_member, # the time step until mask_rim_max. This ensures that for the first t time # steps, the buffer mask keeps increasing. precip_field_mask = ( - precip_forecast_probability_matching_blended >= self.__params.precip_threshold + precip_forecast_probability_matching_blended + >= self.__params.precip_threshold ) # Buffer the mask @@ -2895,7 +3064,9 @@ def __post_process_output(self, j, t_sub, final_blended_forecast_single_member, accumulated_mask = dilated_mask.astype(float) # Iteratively dilate the mask and accumulate the results to create a grayscale rim - mask_rim_temp = min(self.__params.mask_rim + t_sub - 1, self.__params.max_mask_rim) + mask_rim_temp = min( + self.__params.mask_rim + t_sub - 1, self.__params.max_mask_rim + ) for _ in range(mask_rim_temp): dilated_mask = binary_dilation(dilated_mask, struct_element) accumulated_mask += dilated_mask @@ -2905,17 +3076,22 @@ def __post_process_output(self, j, t_sub, final_blended_forecast_single_member, # Get the final mask worker_state.final_blended_forecast_recomposed = ( precip_forecast_min_value - + (worker_state.final_blended_forecast_recomposed - precip_forecast_min_value) + + ( + worker_state.final_blended_forecast_recomposed + - precip_forecast_min_value + ) * precip_field_mask ) precip_field_mask_temp = ( - worker_state.final_blended_forecast_recomposed > precip_forecast_min_value + worker_state.final_blended_forecast_recomposed + > precip_forecast_min_value ) elif self.__config.mask_method == "obs": # The mask equals the most recent benchmark # rainfall field precip_field_mask_temp = ( - precip_forecast_probability_matching_blended >= self.__params.precip_threshold + precip_forecast_probability_matching_blended + >= self.__params.precip_threshold ) # Set to min value outside of mask @@ -2926,18 +3102,23 @@ def __post_process_output(self, j, t_sub, final_blended_forecast_single_member, # If probmatching_method is not None, resample the distribution from # both the extrapolation cascade and the model (NWP) cascade and use # that for the probability matching. - if self.__config.probmatching_method is not None and self.__config.resample_distribution: + if ( + self.__config.probmatching_method is not None + and self.__config.resample_distribution + ): arr1 = worker_state.precip_extrapolated_probability_matching[ worker_state.subtimestep_index ] arr2 = worker_state.precip_models_timestep[j] # resample weights based on cascade level 2. # Areas where one of the fields is nan are not included. - precip_forecast_probability_matching_resampled = probmatching.resample_distributions( - first_array=arr1, - second_array=arr2, - probability_first_array=weights_probability_matching_normalized[0], - randgen=worker_state.randgen_probmatching[j], + precip_forecast_probability_matching_resampled = ( + probmatching.resample_distributions( + first_array=arr1, + second_array=arr2, + probability_first_array=weights_probability_matching_normalized[0], + randgen=worker_state.randgen_probmatching[j], + ) ) else: precip_forecast_probability_matching_resampled = ( @@ -2967,11 +3148,13 @@ def __post_process_output(self, j, t_sub, final_blended_forecast_single_member, # Use R_pm_blended as benchmark field and mean_probabiltity_matching_forecast = np.mean( precip_forecast_probability_matching_resampled[ - precip_forecast_probability_matching_resampled >= self.__params.precip_threshold + precip_forecast_probability_matching_resampled + >= self.__params.precip_threshold ] ) no_rain_mask = ( - worker_state.final_blended_forecast_recomposed >= self.__params.precip_threshold + worker_state.final_blended_forecast_recomposed + >= self.__params.precip_threshold ) mean_precip_forecast = np.mean( worker_state.final_blended_forecast_recomposed[no_rain_mask] @@ -2983,7 +3166,9 @@ def __post_process_output(self, j, t_sub, final_blended_forecast_single_member, ) precip_forecast_probability_matching_resampled = None - final_blended_forecast_single_member.append(worker_state.final_blended_forecast_recomposed) + final_blended_forecast_single_member.append( + worker_state.final_blended_forecast_recomposed + ) return final_blended_forecast_single_member def __measure_time(self, label, start_time): @@ -3533,7 +3718,9 @@ def calculate_weights_spn(correlations, covariance): cov_matrix_inv = inv(cov_matrix) # The component weights are the dot product between cov_matrix_inv and cor_vec weights = np.dot(cov_matrix_inv, correlations) - weights = np.nan_to_num(weights, copy=True, nan=10e-5, posinf=10e-5, neginf=10e-5) + weights = np.nan_to_num( + weights, copy=True, nan=10e-5, posinf=10e-5, neginf=10e-5 + ) weights_dot_correlations = np.dot(weights, correlations) # If the dot product of the weights with the correlations is # larger than 1.0, we assign a weight of 0.0 to the noise (to make @@ -3609,7 +3796,10 @@ def blend_means_sigmas(means, sigmas, weights): sigmas = np.expand_dims(sigmas, axis=sigmas.ndim) # Weight should have one component more (the noise component) than the # means and sigmas. Check this - if weights.shape[0] - means.shape[0] != 1 or weights.shape[0] - sigmas.shape[0] != 1: + if ( + weights.shape[0] - means.shape[0] != 1 + or weights.shape[0] - sigmas.shape[0] != 1 + ): raise ValueError( "The weights array does not have one (noise) component more than mu and sigma" ) From 18a99daa1be0d9b7393fef150908c439bb0e8c98 Mon Sep 17 00:00:00 2001 From: Joep1999 Date: Tue, 23 Sep 2025 09:48:28 +0200 Subject: [PATCH 07/18] fixing bug with probmatching --- pysteps/blending/steps.py | 89 ++++++++++++++++++++++++++++++--------- 1 file changed, 68 insertions(+), 21 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 2013aef8b..d32f168e9 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -797,11 +797,12 @@ def __check_inputs(self): "precip_models must be either a two-dimensional array containing dictionaries with decomposed model fields" "or a four-dimensional array containing the original (NWP) model forecasts" ) - precip_nowcast_dim = self.__precip_nowcast.ndim - if precip_nowcast_dim != 4: - raise ValueError( - "precip_nowcast must be a four-dimensional array containing the externally calculated nowcast" - ) + if self.__precip_nowcast is not None: + precip_nowcast_dim = self.__precip_nowcast.ndim + if precip_nowcast_dim != 4: + raise ValueError( + "precip_nowcast must be a four-dimensional array containing the externally calculated nowcast" + ) if self.__config.extrapolation_kwargs is None: self.__state.extrapolation_kwargs = dict() else: @@ -1045,8 +1046,10 @@ def __prepare_radar_and_NWP_fields(self): # we need to know the zerovalue of precip to replace the mask when decomposing after # extrapolation + self.__params.nowcast_zerovalue = np.nanmin(self.__precip_nowcast) self.__params.precip_zerovalue = np.nanmin(self.__precip) - + print("precip zerovalue is:", self.__params.precip_zerovalue) + print("precip nowcast value is:", self.__params.nowcast_zerovalue) # 1. Start with the radar rainfall fields. We want the fields in a Lagrangian # space. Advect the previous precipitation fields to the same position with # the most recent one (i.e. transform them into the Lagrangian coordinates). @@ -1137,9 +1140,8 @@ def transform_to_lagrangian(precip, i): list(self.__precip_nowcast), ) ) - self.__state.precip_nowcast_cascades = np.array( - [precip_nowcast_decomp] + precip_nowcast_decomp ).swapaxes(1, 2) self.__state.mean_nowcast_timestep = np.array( [precip_nowcast_means] @@ -1752,9 +1754,6 @@ def __find_nowcast_NWP_combination(self, t): With the way it is implemented at this moment: n_ens_members of the output equals the maximum number of (ensemble) members in the input (either the nowcasts or NWP). """ - self.__state.precip_nowcast_timestep = self.__precip_nowcast[:, t, :, :].astype( - np.float64, copy=False - ) self.__state.velocity_models_timestep = self.__velocity_models[ :, t, :, :, : @@ -1766,18 +1765,23 @@ def __find_nowcast_NWP_combination(self, t): "The number of NWP model members is larger than the given number of ensemble members. n_model_members <= n_ens_members." ) - n_ens_members_provided = self.__precip_nowcast.shape[0] - if n_ens_members_provided > self.__config.n_ens_members: - raise ValueError( - "The number of nowcast ensemble members provided is larger than the given number of ensemble members requested. n_ens_members_provided <= n_ens_members." - ) - # Check if NWP models/members should be used individually, or if all of # them are blended together per nowcast ensemble member. if self.__config.blend_nwp_members: self.__state.mapping_list_NWP_member_to_ensemble_member = None elif self.__config.nowcasting_method == "external_nowcast": + + self.__state.precip_nowcast_timestep = self.__precip_nowcast[ + :, t, :, : + ].astype(np.float64, copy=False) + + n_ens_members_provided = self.__precip_nowcast.shape[0] + if n_ens_members_provided > self.__config.n_ens_members: + raise ValueError( + "The number of nowcast ensemble members provided is larger than the given number of ensemble members requested. n_ens_members_provided <= n_ens_members." + ) + n_ens_members_max = max(n_ens_members_provided, n_model_members) n_ens_members_min = min(n_ens_members_provided, n_model_members) print("Number of nowcast ensembles:", n_ens_members_provided) @@ -2258,14 +2262,17 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( if self.__config.nowcasting_method == "external_nowcast": for i in range(self.__config.n_cascade_levels): precip_extrapolated_decomp = worker_state.precip_cascades[j][i][-1, :] + worker_state.time_prev_timestep[j] = t + 1 worker_state.precip_extrapolated_decomp.append( precip_extrapolated_decomp.copy() ) - worker_state.precip_extrapolated_probability_matching.append( - precip_extrapolated_decomp.copy() - ) + + precip_extrapolated = self.__precip_nowcast[j][t][:, :] + worker_state.precip_extrapolated_probability_matching.append( + precip_extrapolated.copy() + ) worker_state.precip_extrapolated_decomp = np.stack( worker_state.precip_extrapolated_decomp @@ -2947,6 +2954,24 @@ def __post_process_output( axis=0, ) # Blend it + rounded_arr = np.round(precip_forecast_probability_matching_final[0], 1) + values, counts = np.unique(rounded_arr, return_counts=True) + + # Find the most common value + most_common = values[np.argmax(counts)] + + print("Most common number before problem var [0]?:", most_common) + print("Count of this number before problem var [0]?:", counts.max()) + + rounded_arr = np.round(precip_forecast_probability_matching_final[1], 1) + values, counts = np.unique(rounded_arr, return_counts=True) + + # Find the most common value + most_common = values[np.argmax(counts)] + + print("Most common number before problem var[1]?:", most_common) + print("Count of this number before problem var[1]?:", counts.max()) + precip_forecast_probability_matching_blended = np.sum( weights_probability_matching_normalized.reshape( weights_probability_matching_normalized.shape[0], 1, 1 @@ -2954,6 +2979,15 @@ def __post_process_output( * precip_forecast_probability_matching_final, axis=0, ) + + rounded_arr = np.round(precip_forecast_probability_matching_blended, 1) + values, counts = np.unique(rounded_arr, return_counts=True) + + # Find the most common value + most_common = values[np.argmax(counts)] + + print("Most common number problem var?:", most_common) + print("Count of this number problem var?:", counts.max()) if self.__config.blend_nwp_members: precip_forecast_probability_matching_blended_mod_only = np.sum( weights_probability_matching_normalized_mod_only.reshape( @@ -3002,7 +3036,16 @@ def __post_process_output( ], axis=0, ) - + print( + "precip_forecast_recomposed_mod_only_no_nan:", + np.nanmin(precip_forecast_recomposed_mod_only_no_nan), + ) + print( + "precip_forecast_recomposed_no_nan:", + np.nanmin(precip_forecast_recomposed_no_nan), + ) + print("mask_radar:", np.nanmin(mask_radar)) + print("mask_model:", np.nanmin(mask_model)) precip_forecast_probability_matching_blended = np.nansum( [ precip_forecast_probability_matching_blended * mask_radar, @@ -3109,6 +3152,7 @@ def __post_process_output( arr1 = worker_state.precip_extrapolated_probability_matching[ worker_state.subtimestep_index ] + print() arr2 = worker_state.precip_models_timestep[j] # resample weights based on cascade level 2. # Areas where one of the fields is nan are not included. @@ -3132,9 +3176,11 @@ def __post_process_output( worker_state.subtimestep_index ] ) + # Adjust the CDF of the forecast to match the resampled distribution combined from # extrapolation and model fields. # Rainfall outside the pure extrapolation domain is not taken into account. + if np.any(np.isfinite(worker_state.final_blended_forecast_recomposed)): worker_state.final_blended_forecast_recomposed = ( probmatching.nonparam_match_empirical_cdf( @@ -3144,6 +3190,7 @@ def __post_process_output( ) ) precip_forecast_probability_matching_resampled = None + elif self.__config.probmatching_method == "mean": # Use R_pm_blended as benchmark field and mean_probabiltity_matching_forecast = np.mean( From f74a60f18e8e6b1d6097e669529a6e22a12c0179 Mon Sep 17 00:00:00 2001 From: Joep1999 Date: Tue, 23 Sep 2025 11:47:34 +0200 Subject: [PATCH 08/18] fixed small problem with means and stds --- pysteps/blending/steps.py | 90 ++++++++++++--------------------------- 1 file changed, 28 insertions(+), 62 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index d32f168e9..515651250 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -1048,8 +1048,6 @@ def __prepare_radar_and_NWP_fields(self): # extrapolation self.__params.nowcast_zerovalue = np.nanmin(self.__precip_nowcast) self.__params.precip_zerovalue = np.nanmin(self.__precip) - print("precip zerovalue is:", self.__params.precip_zerovalue) - print("precip nowcast value is:", self.__params.nowcast_zerovalue) # 1. Start with the radar rainfall fields. We want the fields in a Lagrangian # space. Advect the previous precipitation fields to the same position with # the most recent one (i.e. transform them into the Lagrangian coordinates). @@ -1129,27 +1127,32 @@ def transform_to_lagrangian(precip, i): # Decompose precomputed nowcasts and rearange them again into the required components if self.__precip_nowcast is not None: if self.__precip_nowcast.shape[0] == 1: - precip_nowcast_decomp, precip_nowcast_means, precip_nowcast_stds = ( - self.__decompose_member(self.__precip_nowcast[0]) - ) + results = self.__decompose_member(self.__precip_nowcast[0]) else: with ThreadPool(self.__config.num_workers) as pool: - precip_nowcast_decomp, precip_nowcast_means, precip_nowcast_stds = ( - pool.map( - partial(self.__decompose_member), - list(self.__precip_nowcast), - ) + results = pool.map( + partial(self.__decompose_member), + list(self.__precip_nowcast), ) + precip_nowcast_decomp = [] + precip_nowcast_means = [] + precip_nowcast_stds = [] + + for i in range(self.__precip_nowcast.shape[0]): + precip_nowcast_decomp.append(results[i]["precip_nowcast_decomp"]) + + precip_nowcast_means.append(results[i]["precip_nowcast_means"]) + + precip_nowcast_stds.append(results[i]["precip_nowcast_stds"]) self.__state.precip_nowcast_cascades = np.array( precip_nowcast_decomp ).swapaxes(1, 2) self.__state.mean_nowcast_timestep = np.array( - [precip_nowcast_means] - ).swapaxes(1, 2) - self.__state.std_nowcast_timestep = np.array( - [precip_nowcast_stds] + precip_nowcast_means ).swapaxes(1, 2) - + self.__state.std_nowcast_timestep = np.array(precip_nowcast_stds).swapaxes( + 1, 2 + ) # If necessary, recompose (NWP) model forecasts self.__state.precip_models_cascades = None @@ -1208,8 +1211,13 @@ def __decompose_member(self, member_field): results_decomp.append(res["cascade_levels"]) means.append(res["means"]) stds.append(res["stds"]) + results = { + "precip_nowcast_decomp": results_decomp, + "precip_nowcast_means": means, + "precip_nowcast_stds": stds, + } - return results_decomp, means, stds + return results def __zero_precipitation_forecast(self): """ @@ -1784,8 +1792,6 @@ def __find_nowcast_NWP_combination(self, t): n_ens_members_max = max(n_ens_members_provided, n_model_members) n_ens_members_min = min(n_ens_members_provided, n_model_members) - print("Number of nowcast ensembles:", n_ens_members_provided) - print("Number of nwp ensembles:", n_model_members) # Also make a list of the model index numbers. These indices are needed # for indexing the right climatological skill file when pysteps calculates # the blended forecast in parallel. @@ -2186,7 +2192,6 @@ def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): # If nowcast method seleced is external_nowcast, n_ens members has to be 1 if self.__config.nowcasting_method == "external_nowcast": - print("Using nowcasting method:", self.__config.nowcasting_method) for i in range(self.__config.n_cascade_levels): # Use a deterministic Externally computed nowcasting model worker_state.precip_cascades[j][i] = ( @@ -2195,7 +2200,6 @@ def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): # Follow the 'standard' STEPS blending approach as described in :cite:`Imhoff2023` elif self.__config.nowcasting_method == "steps": - print("Using nowcasting method:", self.__config.nowcasting_method) for i in range(self.__config.n_cascade_levels): # apply AR(p) process to extrapolation cascade level if ( @@ -2659,14 +2663,14 @@ def __blend_cascades(self, t_sub, j, worker_state): ) # [(extr_field, n_model_fields), n_cascade_levels, ...] means_stacked = np.concatenate( ( - worker_state.mean_nowcast_timestep[t_sub, j], + worker_state.mean_nowcast_timestep[None, j, :, t_sub], worker_state.mean_models_timestep, ), axis=0, ) sigmas_stacked = np.concatenate( ( - worker_state.std_nowcast_timestep[t_sub, j], + worker_state.std_nowcast_timestep[None, j, :, t_sub], worker_state.std_models_timestep, ), axis=0, @@ -2716,7 +2720,7 @@ def __blend_cascades(self, t_sub, j, worker_state): means_stacked = np.concatenate( ( # worker_state.mean_nowcast_timestep, - worker_state.mean_nowcast_timestep[t_sub, j], + worker_state.mean_nowcast_timestep[None, j, :, t_sub], worker_state.mean_models_timestep[None, j], ), axis=0, @@ -2724,7 +2728,7 @@ def __blend_cascades(self, t_sub, j, worker_state): sigmas_stacked = np.concatenate( ( # worker_state.std_nowcast_timestep, - worker_state.std_nowcast_timestep[t_sub, j], + worker_state.std_nowcast_timestep[None, j, :, t_sub], worker_state.std_models_timestep[None, j], ), axis=0, @@ -2954,24 +2958,6 @@ def __post_process_output( axis=0, ) # Blend it - rounded_arr = np.round(precip_forecast_probability_matching_final[0], 1) - values, counts = np.unique(rounded_arr, return_counts=True) - - # Find the most common value - most_common = values[np.argmax(counts)] - - print("Most common number before problem var [0]?:", most_common) - print("Count of this number before problem var [0]?:", counts.max()) - - rounded_arr = np.round(precip_forecast_probability_matching_final[1], 1) - values, counts = np.unique(rounded_arr, return_counts=True) - - # Find the most common value - most_common = values[np.argmax(counts)] - - print("Most common number before problem var[1]?:", most_common) - print("Count of this number before problem var[1]?:", counts.max()) - precip_forecast_probability_matching_blended = np.sum( weights_probability_matching_normalized.reshape( weights_probability_matching_normalized.shape[0], 1, 1 @@ -2979,15 +2965,6 @@ def __post_process_output( * precip_forecast_probability_matching_final, axis=0, ) - - rounded_arr = np.round(precip_forecast_probability_matching_blended, 1) - values, counts = np.unique(rounded_arr, return_counts=True) - - # Find the most common value - most_common = values[np.argmax(counts)] - - print("Most common number problem var?:", most_common) - print("Count of this number problem var?:", counts.max()) if self.__config.blend_nwp_members: precip_forecast_probability_matching_blended_mod_only = np.sum( weights_probability_matching_normalized_mod_only.reshape( @@ -3036,16 +3013,6 @@ def __post_process_output( ], axis=0, ) - print( - "precip_forecast_recomposed_mod_only_no_nan:", - np.nanmin(precip_forecast_recomposed_mod_only_no_nan), - ) - print( - "precip_forecast_recomposed_no_nan:", - np.nanmin(precip_forecast_recomposed_no_nan), - ) - print("mask_radar:", np.nanmin(mask_radar)) - print("mask_model:", np.nanmin(mask_model)) precip_forecast_probability_matching_blended = np.nansum( [ precip_forecast_probability_matching_blended * mask_radar, @@ -3152,7 +3119,6 @@ def __post_process_output( arr1 = worker_state.precip_extrapolated_probability_matching[ worker_state.subtimestep_index ] - print() arr2 = worker_state.precip_models_timestep[j] # resample weights based on cascade level 2. # Areas where one of the fields is nan are not included. From e6beb32ee16f8f810116b8d68a2120d59a7033f5 Mon Sep 17 00:00:00 2001 From: Ruben Imhoff Date: Tue, 23 Sep 2025 11:52:41 +0200 Subject: [PATCH 09/18] feat: first steps towards allowing noise development with external nowcast --- pysteps/blending/steps.py | 403 ++++++++++++++++++++------------------ 1 file changed, 208 insertions(+), 195 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 2013aef8b..7aafc01ca 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -1090,11 +1090,13 @@ def transform_to_lagrangian(precip, i): if self.__precip_nowcast is not None: self.__precip_nowcast = self.__precip_nowcast.copy() # TODO: check this function for ensemble nowcasts -> now made it always work with 4 dimentions - for m in range(self.__precip_nowcast.shape[0]): + for ens_mem in range(self.__precip_nowcast.shape[0]): for t in range(self.__precip_nowcast.shape[1]): self.__precip_nowcast[ - m, t, ~np.isfinite(self.__precip_nowcast[m, t, :, :]) - ] = np.nanmin(self.__precip_nowcast[m, t, :, :]) + ens_mem, + t, + ~np.isfinite(self.__precip_nowcast[ens_mem, t, :, :]), + ] = np.nanmin(self.__precip_nowcast[ens_mem, t, :, :]) # Perform the cascade decomposition for the input precip fields and, # if necessary, for the (NWP) model fields @@ -1175,7 +1177,6 @@ def transform_to_lagrangian(precip, i): self.__params.noise_kwargs["win_fun"], ) - # TODO This variable is currently not used anywhere. Should we calculate it? # The norain fraction threshold used for nwp is the default value of 0.0, # since nwp does not suffer from clutter. self.__params.zero_precip_model_fields = check_norain( @@ -2179,8 +2180,6 @@ def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): # Regress the extrapolation component to the subsequent time step. # Iterate the AR(p) model for each cascade level - - # If nowcast method seleced is external_nowcast, n_ens members has to be 1 if self.__config.nowcasting_method == "external_nowcast": print("Using nowcasting method:", self.__config.nowcasting_method) for i in range(self.__config.n_cascade_levels): @@ -2275,12 +2274,12 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( worker_state.precip_extrapolated_probability_matching ) # [None, :] - else: - # Extrapolate per sub time step - for t_sub in worker_state.subtimesteps: - if t_sub > 0: - t_diff_prev_subtimestep_int = t_sub - int(t_sub) - if t_diff_prev_subtimestep_int > 0.0: + # Extrapolate per sub time step + for t_sub in worker_state.subtimesteps: + if t_sub > 0: + t_diff_prev_subtimestep_int = t_sub - int(t_sub) + if t_diff_prev_subtimestep_int > 0.0: + if self.__config.nowcasting_method == "steps": precip_forecast_cascade_subtimestep = [ (1.0 - t_diff_prev_subtimestep_int) * worker_state.precip_cascades_prev_subtimestep[j][i][-1, :] @@ -2288,6 +2287,7 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( * worker_state.precip_cascades[j][i][-1, :] for i in range(self.__config.n_cascade_levels) ] + if self.__config.noise_method is not None: noise_cascade_subtimestep = [ (1.0 - t_diff_prev_subtimestep_int) * worker_state.cascade_noise_prev_subtimestep[j][i][-1, :] @@ -2296,74 +2296,79 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( for i in range(self.__config.n_cascade_levels) ] - else: + else: + if self.__config.nowcasting_method == "steps": precip_forecast_cascade_subtimestep = [ worker_state.precip_cascades_prev_subtimestep[j][i][-1, :] for i in range(self.__config.n_cascade_levels) ] + if self.__config.noise_method is not None: noise_cascade_subtimestep = [ worker_state.cascade_noise_prev_subtimestep[j][i][-1, :] for i in range(self.__config.n_cascade_levels) ] + if self.__config.nowcasting_method == "steps": precip_forecast_cascade_subtimestep = np.stack( precip_forecast_cascade_subtimestep ) + if self.__config.noise_method is not None: noise_cascade_subtimestep = np.stack(noise_cascade_subtimestep) - t_diff_prev_subtimestep = t_sub - worker_state.time_prev_timestep[j] - worker_state.leadtime_since_start_forecast[ - j - ] += t_diff_prev_subtimestep - - # compute the perturbed motion field - include the NWP - # velocities and the weights. Note that we only perturb - # the extrapolation velocity field, as the NWP velocity - # field is present per time step - if self.__config.velocity_perturbation_method is not None: - velocity_perturbations_extrapolation = ( - self.__velocity - + self.__params.generate_velocity_noise( - self.__params.velocity_perturbations[j], - worker_state.leadtime_since_start_forecast[j] - * self.__config.timestep, - ) - ) + t_diff_prev_subtimestep = t_sub - worker_state.time_prev_timestep[j] + worker_state.leadtime_since_start_forecast[j] += t_diff_prev_subtimestep - # Stack the perturbed extrapolation and the NWP velocities - if self.__config.blend_nwp_members: - velocity_stack_all = np.concatenate( - ( - velocity_perturbations_extrapolation[None, :, :, :], - worker_state.velocity_models_timestep, - ), - axis=0, - ) - else: - velocity_models = worker_state.velocity_models_timestep[j] - velocity_stack_all = np.concatenate( - ( - velocity_perturbations_extrapolation[None, :, :, :], - velocity_models[None, :, :, :], - ), - axis=0, + # compute the perturbed motion field - include the NWP + # velocities and the weights. Note that we only perturb + # the extrapolation velocity field, as the NWP velocity + # field is present per time step + if self.__config.velocity_perturbation_method is not None: + velocity_perturbations_extrapolation = ( + self.__velocity + + self.__params.generate_velocity_noise( + self.__params.velocity_perturbations[j], + worker_state.leadtime_since_start_forecast[j] + * self.__config.timestep, ) - velocity_models = None - - # Obtain a blended optical flow, using the weights of the - # second cascade following eq. 24 in BPS2006 - velocity_blended = blending.utils.blend_optical_flows( - flows=velocity_stack_all, - weights=worker_state.weights[ - :-1, 1 - ], # [(extr_field, n_model_fields), cascade_level=2] - ) - - # Extrapolate both cascades to the next time step - # First recompose the cascade, advect it and decompose it again - # This is needed to remove the interpolation artifacts. - # In addition, the number of extrapolations is greatly reduced - # A. Radar Rain + ) + + # Stack the perturbed extrapolation and the NWP velocities + if self.__config.blend_nwp_members: + velocity_stack_all = np.concatenate( + ( + velocity_perturbations_extrapolation[None, :, :, :], + worker_state.velocity_models_timestep, + ), + axis=0, + ) + else: + velocity_models = worker_state.velocity_models_timestep[j] + velocity_stack_all = np.concatenate( + ( + velocity_perturbations_extrapolation[None, :, :, :], + velocity_models[None, :, :, :], + ), + axis=0, + ) + velocity_models = None + + # Obtain a blended optical flow, using the weights of the + # second cascade following eq. 24 in BPS2006 + velocity_blended = blending.utils.blend_optical_flows( + flows=velocity_stack_all, + weights=worker_state.weights[ + :-1, 1 + ], # [(extr_field, n_model_fields), cascade_level=2] + ) + + # Extrapolate both cascades to the next time step + # First recompose the cascade, advect it and decompose it again + # This is needed to remove the interpolation artifacts. + # In addition, the number of extrapolations is greatly reduced + + # A. The extrapolation component + if self.__config.nowcasting_method == "steps": + # First, recompose the cascades into one forecast precip_forecast_recomp_subtimestep = ( blending.utils.recompose_cascade( combined_cascade=precip_forecast_cascade_subtimestep, @@ -2401,10 +2406,11 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( precip_forecast_extrapolated_recomp_subtimestep_temp[0].copy() ) temp_mask = ~np.isfinite(precip_extrapolated_recomp_subtimestep) - # TODO: WHERE DO CAN I FIND THIS -15.0 + # Set non-finite values to the zerovalue precip_extrapolated_recomp_subtimestep[ ~np.isfinite(precip_extrapolated_recomp_subtimestep) ] = self.__params.precip_zerovalue + # Decompose the forecast again into multiplicative cascades precip_extrapolated_decomp = self.__params.decomposition_method( precip_extrapolated_recomp_subtimestep, self.__params.bandpass_filter, @@ -2426,7 +2432,21 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( ) for i in range(self.__config.n_cascade_levels): precip_extrapolated_decomp[i][temp_mask] = np.nan - # B. Noise + + # Append the results to the output lists + worker_state.precip_extrapolated_decomp.append( + precip_extrapolated_decomp.copy() + ) + + precip_forecast_cascade_subtimestep = None + precip_forecast_recomp_subtimestep = None + precip_forecast_extrapolated_recomp_subtimestep_temp = None + precip_extrapolated_recomp_subtimestep = None + precip_extrapolated_decomp = None + + # B. The noise component + if self.__config.noise_method is not None: + # First, recompose the cascades into one forecast noise_cascade_subtimestep_recomp = blending.utils.recompose_cascade( combined_cascade=noise_cascade_subtimestep, combined_mean=worker_state.precip_mean_noise[j], @@ -2447,6 +2467,7 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( **extrap_kwargs_noise, ) noise_extrapolated_recomp = noise_extrapolated_recomp_temp[0].copy() + # Decompose the noise component again into multiplicative cascades noise_extrapolated_decomp = self.__params.decomposition_method( noise_extrapolated_recomp, self.__params.bandpass_filter, @@ -2463,159 +2484,151 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( ] # Append the results to the output lists - worker_state.precip_extrapolated_decomp.append( - precip_extrapolated_decomp.copy() - ) worker_state.noise_extrapolated_decomp.append( noise_extrapolated_decomp.copy() ) - precip_forecast_cascade_subtimestep = None - precip_forecast_recomp_subtimestep = None - precip_forecast_extrapolated_recomp_subtimestep_temp = None - precip_extrapolated_recomp_subtimestep = None - precip_extrapolated_decomp = None + noise_cascade_subtimestep = None noise_cascade_subtimestep_recomp = None noise_extrapolated_recomp_temp = None noise_extrapolated_recomp = None noise_extrapolated_decomp = None - # Finally, also extrapolate the initial radar rainfall - # field. This will be blended with the rainfall field(s) - # of the (NWP) model(s) for Lagrangian blended prob. matching - # min_R = np.min(precip) - extrap_kwargs_pb["displacement_prev"] = ( - worker_state.previous_displacement_prob_matching[j] - ) - # Apply the domain mask to the extrapolation component - precip_forecast_temp_for_probability_matching = self.__precip.copy() - precip_forecast_temp_for_probability_matching[ - self.__params.domain_mask - ] = np.nan - - ( - precip_forecast_extrapolated_probability_matching_temp, - worker_state.previous_displacement_prob_matching[j], - ) = self.__params.extrapolation_method( - precip_forecast_temp_for_probability_matching, - velocity_blended, - [t_diff_prev_subtimestep], - allow_nonfinite_values=True, - **extrap_kwargs_pb, - ) - - worker_state.precip_extrapolated_probability_matching.append( - precip_forecast_extrapolated_probability_matching_temp[0] - ) - - worker_state.time_prev_timestep[j] = t_sub - - if len(worker_state.precip_extrapolated_decomp) > 0: - worker_state.precip_extrapolated_decomp = np.stack( - worker_state.precip_extrapolated_decomp + # Finally, also extrapolate the initial radar rainfall + # field. This will be blended with the rainfall field(s) + # of the (NWP) model(s) for Lagrangian blended prob. matching + # min_R = np.min(precip) + extrap_kwargs_pb["displacement_prev"] = ( + worker_state.previous_displacement_prob_matching[j] ) - worker_state.noise_extrapolated_decomp = np.stack( - worker_state.noise_extrapolated_decomp + # Apply the domain mask to the extrapolation component + precip_forecast_temp_for_probability_matching = self.__precip.copy() + precip_forecast_temp_for_probability_matching[ + self.__params.domain_mask + ] = np.nan + + ( + precip_forecast_extrapolated_probability_matching_temp, + worker_state.previous_displacement_prob_matching[j], + ) = self.__params.extrapolation_method( + precip_forecast_temp_for_probability_matching, + velocity_blended, + [t_diff_prev_subtimestep], + allow_nonfinite_values=True, + **extrap_kwargs_pb, ) - worker_state.precip_extrapolated_probability_matching = np.stack( - worker_state.precip_extrapolated_probability_matching + + worker_state.precip_extrapolated_probability_matching.append( + precip_forecast_extrapolated_probability_matching_temp[0] ) - # advect the forecast field by one time step if no subtimesteps in the - # current interval were found - if not worker_state.subtimesteps: - t_diff_prev_subtimestep = t + 1 - worker_state.time_prev_timestep[j] - worker_state.leadtime_since_start_forecast[j] += t_diff_prev_subtimestep + worker_state.time_prev_timestep[j] = t_sub + # TODO: continue here - # compute the perturbed motion field - include the NWP - # velocities and the weights - if self.__config.velocity_perturbation_method is not None: - velocity_perturbations_extrapolation = ( - self.__velocity - + self.__params.generate_velocity_noise( - self.__params.velocity_perturbations[j], - worker_state.leadtime_since_start_forecast[j] - * self.__config.timestep, - ) - ) + if len(worker_state.precip_extrapolated_decomp) > 0: + worker_state.precip_extrapolated_decomp = np.stack( + worker_state.precip_extrapolated_decomp + ) + worker_state.noise_extrapolated_decomp = np.stack( + worker_state.noise_extrapolated_decomp + ) + worker_state.precip_extrapolated_probability_matching = np.stack( + worker_state.precip_extrapolated_probability_matching + ) - # Stack the perturbed extrapolation and the NWP velocities - if self.__config.blend_nwp_members: - velocity_stack_all = np.concatenate( - ( - velocity_perturbations_extrapolation[None, :, :, :], - worker_state.velocity_models_timestep, - ), - axis=0, - ) - else: - velocity_models = worker_state.velocity_models_timestep[j] - velocity_stack_all = np.concatenate( - ( - velocity_perturbations_extrapolation[None, :, :, :], - velocity_models[None, :, :, :], - ), - axis=0, - ) - velocity_models = None + # advect the forecast field by one time step if no subtimesteps in the + # current interval were found + if not worker_state.subtimesteps: + t_diff_prev_subtimestep = t + 1 - worker_state.time_prev_timestep[j] + worker_state.leadtime_since_start_forecast[j] += t_diff_prev_subtimestep - # Obtain a blended optical flow, using the weights of the - # second cascade following eq. 24 in BPS2006 - velocity_blended = blending.utils.blend_optical_flows( - flows=velocity_stack_all, - weights=worker_state.weights[ - :-1, 1 - ], # [(extr_field, n_model_fields), cascade_level=2] + # compute the perturbed motion field - include the NWP + # velocities and the weights + if self.__config.velocity_perturbation_method is not None: + velocity_perturbations_extrapolation = ( + self.__velocity + + self.__params.generate_velocity_noise( + self.__params.velocity_perturbations[j], + worker_state.leadtime_since_start_forecast[j] + * self.__config.timestep, + ) ) - # Extrapolate the extrapolation and noise cascade - extrap_kwargs_["displacement_prev"] = ( - worker_state.previous_displacement[j] + # Stack the perturbed extrapolation and the NWP velocities + if self.__config.blend_nwp_members: + velocity_stack_all = np.concatenate( + ( + velocity_perturbations_extrapolation[None, :, :, :], + worker_state.velocity_models_timestep, + ), + axis=0, ) - extrap_kwargs_noise["displacement_prev"] = ( - worker_state.previous_displacement_noise_cascade[j] + else: + velocity_models = worker_state.velocity_models_timestep[j] + velocity_stack_all = np.concatenate( + ( + velocity_perturbations_extrapolation[None, :, :, :], + velocity_models[None, :, :, :], + ), + axis=0, ) - extrap_kwargs_noise["map_coordinates_mode"] = "wrap" + velocity_models = None - ( - _, - worker_state.previous_displacement[j], - ) = self.__params.extrapolation_method( - None, - velocity_blended, - [t_diff_prev_subtimestep], - allow_nonfinite_values=True, - **extrap_kwargs_, - ) + # Obtain a blended optical flow, using the weights of the + # second cascade following eq. 24 in BPS2006 + velocity_blended = blending.utils.blend_optical_flows( + flows=velocity_stack_all, + weights=worker_state.weights[ + :-1, 1 + ], # [(extr_field, n_model_fields), cascade_level=2] + ) - ( - _, - worker_state.previous_displacement_noise_cascade[j], - ) = self.__params.extrapolation_method( - None, - velocity_blended, - [t_diff_prev_subtimestep], - allow_nonfinite_values=True, - **extrap_kwargs_noise, - ) + # Extrapolate the extrapolation and noise cascade + extrap_kwargs_["displacement_prev"] = worker_state.previous_displacement[j] + extrap_kwargs_noise["displacement_prev"] = ( + worker_state.previous_displacement_noise_cascade[j] + ) + extrap_kwargs_noise["map_coordinates_mode"] = "wrap" - # Also extrapolate the radar observation, used for the probability - # matching and post-processing steps - extrap_kwargs_pb["displacement_prev"] = ( - worker_state.previous_displacement_prob_matching[j] - ) - ( - _, - worker_state.previous_displacement_prob_matching[j], - ) = self.__params.extrapolation_method( - None, - velocity_blended, - [t_diff_prev_subtimestep], - allow_nonfinite_values=True, - **extrap_kwargs_pb, - ) + ( + _, + worker_state.previous_displacement[j], + ) = self.__params.extrapolation_method( + None, + velocity_blended, + [t_diff_prev_subtimestep], + allow_nonfinite_values=True, + **extrap_kwargs_, + ) - worker_state.time_prev_timestep[j] = t + 1 + ( + _, + worker_state.previous_displacement_noise_cascade[j], + ) = self.__params.extrapolation_method( + None, + velocity_blended, + [t_diff_prev_subtimestep], + allow_nonfinite_values=True, + **extrap_kwargs_noise, + ) + + # Also extrapolate the radar observation, used for the probability + # matching and post-processing steps + extrap_kwargs_pb["displacement_prev"] = ( + worker_state.previous_displacement_prob_matching[j] + ) + ( + _, + worker_state.previous_displacement_prob_matching[j], + ) = self.__params.extrapolation_method( + None, + velocity_blended, + [t_diff_prev_subtimestep], + allow_nonfinite_values=True, + **extrap_kwargs_pb, + ) + + worker_state.time_prev_timestep[j] = t + 1 worker_state.precip_cascades_prev_subtimestep[j] = worker_state.precip_cascades[ j From ff7108495b468bbe8720d8793393b033e7a1df53 Mon Sep 17 00:00:00 2001 From: Ruben Imhoff Date: Tue, 23 Sep 2025 13:23:58 +0200 Subject: [PATCH 10/18] feat: allow adding noise to external nowcast blending --- pysteps/blending/steps.py | 272 ++++++++++++++++++++++---------------- 1 file changed, 155 insertions(+), 117 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 7198c92b8..f725d729b 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -1780,7 +1780,6 @@ def __find_nowcast_NWP_combination(self, t): self.__state.mapping_list_NWP_member_to_ensemble_member = None elif self.__config.nowcasting_method == "external_nowcast": - self.__state.precip_nowcast_timestep = self.__precip_nowcast[ :, t, :, : ].astype(np.float64, copy=False) @@ -2262,29 +2261,6 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( worker_state.noise_extrapolated_decomp = [] worker_state.precip_extrapolated_probability_matching = [] - if self.__config.nowcasting_method == "external_nowcast": - for i in range(self.__config.n_cascade_levels): - precip_extrapolated_decomp = worker_state.precip_cascades[j][i][-1, :] - - worker_state.time_prev_timestep[j] = t + 1 - - worker_state.precip_extrapolated_decomp.append( - precip_extrapolated_decomp.copy() - ) - - precip_extrapolated = self.__precip_nowcast[j][t][:, :] - worker_state.precip_extrapolated_probability_matching.append( - precip_extrapolated.copy() - ) - - worker_state.precip_extrapolated_decomp = np.stack( - worker_state.precip_extrapolated_decomp - )[None, :] - - worker_state.precip_extrapolated_probability_matching = np.stack( - worker_state.precip_extrapolated_probability_matching - ) # [None, :] - # Extrapolate per sub time step for t_sub in worker_state.subtimesteps: if t_sub > 0: @@ -2374,7 +2350,7 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( # Extrapolate both cascades to the next time step # First recompose the cascade, advect it and decompose it again - # This is needed to remove the interpolation artifacts. + # This is needed to remove the interpolation artefacts. # In addition, the number of extrapolations is greatly reduced # A. The extrapolation component @@ -2505,47 +2481,49 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( noise_extrapolated_recomp = None noise_extrapolated_decomp = None - # Finally, also extrapolate the initial radar rainfall - # field. This will be blended with the rainfall field(s) - # of the (NWP) model(s) for Lagrangian blended prob. matching - # min_R = np.min(precip) - extrap_kwargs_pb["displacement_prev"] = ( - worker_state.previous_displacement_prob_matching[j] - ) - # Apply the domain mask to the extrapolation component - precip_forecast_temp_for_probability_matching = self.__precip.copy() - precip_forecast_temp_for_probability_matching[ - self.__params.domain_mask - ] = np.nan + # Finally, also extrapolate the initial radar rainfall field. This will be + # blended with the rainfall field(s) of the (NWP) model(s) for Lagrangian + # blended prob. matching min_R = np.min(precip). If we use an external + # nowcast, this variable will be set later in this function. + if self.__config.nowcasting_method == "steps": + extrap_kwargs_pb["displacement_prev"] = ( + worker_state.previous_displacement_prob_matching[j] + ) + # Apply the domain mask to the extrapolation component + precip_forecast_temp_for_probability_matching = self.__precip.copy() + precip_forecast_temp_for_probability_matching[ + self.__params.domain_mask + ] = np.nan - ( - precip_forecast_extrapolated_probability_matching_temp, - worker_state.previous_displacement_prob_matching[j], - ) = self.__params.extrapolation_method( - precip_forecast_temp_for_probability_matching, - velocity_blended, - [t_diff_prev_subtimestep], - allow_nonfinite_values=True, - **extrap_kwargs_pb, - ) + ( + precip_forecast_extrapolated_probability_matching_temp, + worker_state.previous_displacement_prob_matching[j], + ) = self.__params.extrapolation_method( + precip_forecast_temp_for_probability_matching, + velocity_blended, + [t_diff_prev_subtimestep], + allow_nonfinite_values=True, + **extrap_kwargs_pb, + ) - worker_state.precip_extrapolated_probability_matching.append( - precip_forecast_extrapolated_probability_matching_temp[0] - ) + worker_state.precip_extrapolated_probability_matching.append( + precip_forecast_extrapolated_probability_matching_temp[0] + ) - worker_state.time_prev_timestep[j] = t_sub - # TODO: continue here + worker_state.time_prev_timestep[j] = t_sub if len(worker_state.precip_extrapolated_decomp) > 0: - worker_state.precip_extrapolated_decomp = np.stack( - worker_state.precip_extrapolated_decomp - ) - worker_state.noise_extrapolated_decomp = np.stack( - worker_state.noise_extrapolated_decomp - ) - worker_state.precip_extrapolated_probability_matching = np.stack( - worker_state.precip_extrapolated_probability_matching - ) + if self.__config.nowcasting_method == "steps": + worker_state.precip_extrapolated_decomp = np.stack( + worker_state.precip_extrapolated_decomp + ) + worker_state.precip_extrapolated_probability_matching = np.stack( + worker_state.precip_extrapolated_probability_matching + ) + if self.__config.noise_method is not None: + worker_state.noise_extrapolated_decomp = np.stack( + worker_state.noise_extrapolated_decomp + ) # advect the forecast field by one time step if no subtimesteps in the # current interval were found @@ -2601,46 +2579,77 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( ) extrap_kwargs_noise["map_coordinates_mode"] = "wrap" - ( - _, - worker_state.previous_displacement[j], - ) = self.__params.extrapolation_method( - None, - velocity_blended, - [t_diff_prev_subtimestep], - allow_nonfinite_values=True, - **extrap_kwargs_, - ) - - ( - _, - worker_state.previous_displacement_noise_cascade[j], - ) = self.__params.extrapolation_method( - None, - velocity_blended, - [t_diff_prev_subtimestep], - allow_nonfinite_values=True, - **extrap_kwargs_noise, - ) - + # Extrapolate the extrapolation cascade + if self.__config.nowcasting_method == "steps": + ( + _, + worker_state.previous_displacement[j], + ) = self.__params.extrapolation_method( + None, + velocity_blended, + [t_diff_prev_subtimestep], + allow_nonfinite_values=True, + **extrap_kwargs_, + ) + # Extrapolate the noise cascade + if self.__config.noise_method is not None: + ( + _, + worker_state.previous_displacement_noise_cascade[j], + ) = self.__params.extrapolation_method( + None, + velocity_blended, + [t_diff_prev_subtimestep], + allow_nonfinite_values=True, + **extrap_kwargs_noise, + ) # Also extrapolate the radar observation, used for the probability # matching and post-processing steps - extrap_kwargs_pb["displacement_prev"] = ( - worker_state.previous_displacement_prob_matching[j] - ) - ( - _, - worker_state.previous_displacement_prob_matching[j], - ) = self.__params.extrapolation_method( - None, - velocity_blended, - [t_diff_prev_subtimestep], - allow_nonfinite_values=True, - **extrap_kwargs_pb, - ) + if self.__config.nowcasting_method == "steps": + extrap_kwargs_pb["displacement_prev"] = ( + worker_state.previous_displacement_prob_matching[j] + ) + ( + _, + worker_state.previous_displacement_prob_matching[j], + ) = self.__params.extrapolation_method( + None, + velocity_blended, + [t_diff_prev_subtimestep], + allow_nonfinite_values=True, + **extrap_kwargs_pb, + ) worker_state.time_prev_timestep[j] = t + 1 + # If an external nowcast is provided, precip_extrapolated_decomp and + # precip_extrapolated_probability_matching have been omitted so far. + # Fill them in with the external nowcast information now. + if self.__config.nowcasting_method == "external_nowcast": + for i in range(self.__config.n_cascade_levels): + precip_extrapolated_decomp = worker_state.precip_cascades[j][i][-1, :] + + worker_state.time_prev_timestep[j] = t + 1 + + worker_state.precip_extrapolated_decomp.append( + precip_extrapolated_decomp.copy() + ) + + # Also update the probability matching fields + precip_extrapolated = self.__precip_nowcast[j][t][:, :] + worker_state.precip_extrapolated_probability_matching.append( + precip_extrapolated.copy() + ) + + # Stack it for the output + worker_state.precip_extrapolated_decomp = np.stack( + worker_state.precip_extrapolated_decomp + )[None, :] + + worker_state.precip_extrapolated_probability_matching = np.stack( + worker_state.precip_extrapolated_probability_matching + ) # [None, :] + worker_state.precip_cascades_prev_subtimestep[j] = worker_state.precip_cascades[ j ] @@ -2659,21 +2668,33 @@ def __blend_cascades(self, t_sub, j, worker_state): )[0][0] # First concatenate the cascades and the means and sigmas # precip_models = [n_models,timesteps,n_cascade_levels,m,n] - - # If external_nowcast is used, do not use and blend with the noise component if ( self.__config.blend_nwp_members and self.__config.nowcasting_method == "external_nowcast" ): - cascade_stack_all_components = np.concatenate( - ( - worker_state.precip_extrapolated_decomp[ - None, worker_state.subtimestep_index - ], - worker_state.precip_models_cascades_timestep, - ), - axis=0, - ) # [(extr_field, n_model_fields), n_cascade_levels, ...] + if self.__config.noise_method is None: + cascade_stack_all_components = np.concatenate( + ( + worker_state.precip_extrapolated_decomp[ + None, worker_state.subtimestep_index + ], + worker_state.precip_models_cascades_timestep, + ), + axis=0, + ) # [(extr_field, n_model_fields), n_cascade_levels, ...] + else: + cascade_stack_all_components = np.concatenate( + ( + worker_state.precip_extrapolated_decomp[ + None, worker_state.subtimestep_index + ], + worker_state.precip_models_cascades_timestep, + worker_state.noise_extrapolated_decomp[ + None, worker_state.subtimestep_index + ], + ), + axis=0, + ) # [(extr_field, n_model_fields), n_cascade_levels, ...] means_stacked = np.concatenate( ( worker_state.mean_nowcast_timestep[None, j, :, t_sub], @@ -2721,15 +2742,29 @@ def __blend_cascades(self, t_sub, j, worker_state): ) elif self.__config.nowcasting_method == "external_nowcast": - cascade_stack_all_components = np.concatenate( - ( - worker_state.precip_extrapolated_decomp[ - None, worker_state.subtimestep_index - ], - worker_state.precip_models_cascades_timestep[None, j], - ), - axis=0, - ) # [(extr_field, n_model_fields), n_cascade_levels, ...] + if self.__config.noise_method is None: + cascade_stack_all_components = np.concatenate( + ( + worker_state.precip_extrapolated_decomp[ + None, worker_state.subtimestep_index + ], + worker_state.precip_models_cascades_timestep[None, j], + ), + axis=0, + ) # [(extr_field, n_model_fields), n_cascade_levels, ...] + else: + cascade_stack_all_components = np.concatenate( + ( + worker_state.precip_extrapolated_decomp[ + None, worker_state.subtimestep_index + ], + worker_state.precip_models_cascades_timestep[None, j], + worker_state.noise_extrapolated_decomp[ + None, worker_state.subtimestep_index + ], + ), + axis=0, + ) # [(extr_field, n_model_fields), n_cascade_levels, ...] means_stacked = np.concatenate( ( # worker_state.mean_nowcast_timestep, @@ -2808,7 +2843,10 @@ def __blend_cascades(self, t_sub, j, worker_state): worker_state.weights_model_only_with_noise = ( worker_state.weights_model_only.copy() ) - if self.__config.nowcasting_method == "external_nowcast": + if ( + self.__config.nowcasting_method == "external_nowcast" + and self.__config.noise_method is None + ): # First determine the weights without noise worker_state.weights = worker_state.weights[:-1, :] / np.sum( worker_state.weights[:-1, :], axis=0 From 31628c02408c58b99468ac19a190d2301142425e Mon Sep 17 00:00:00 2001 From: Joep1999 Date: Tue, 23 Sep 2025 13:46:28 +0200 Subject: [PATCH 11/18] added support for copying multiple nowcasts and forecasts in case n_ens_asked is bigger --- pysteps/blending/steps.py | 242 ++++++++++++++++++++++++-------------- 1 file changed, 151 insertions(+), 91 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index f725d729b..9a218f935 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -1130,31 +1130,44 @@ def transform_to_lagrangian(precip, i): if self.__precip_nowcast is not None: if self.__precip_nowcast.shape[0] == 1: results = self.__decompose_member(self.__precip_nowcast[0]) + self.__state.precip_nowcast_cascades = np.array( + results["precip_nowcast_decomp"] + ).swapaxes(0, 1) + self.__state.mean_nowcast_timestep = np.array( + results["precip_nowcast_means"] + ).swapaxes(0, 1) + self.__state.std_nowcast_timestep = np.array( + results["precip_nowcast_stds"] + ).swapaxes(0, 1) + print( + "precip_nowcast_cascades shape:", + self.__state.precip_nowcast_cascades.shape, + ) else: with ThreadPool(self.__config.num_workers) as pool: results = pool.map( partial(self.__decompose_member), list(self.__precip_nowcast), ) - precip_nowcast_decomp = [] - precip_nowcast_means = [] - precip_nowcast_stds = [] - - for i in range(self.__precip_nowcast.shape[0]): - precip_nowcast_decomp.append(results[i]["precip_nowcast_decomp"]) - - precip_nowcast_means.append(results[i]["precip_nowcast_means"]) - - precip_nowcast_stds.append(results[i]["precip_nowcast_stds"]) - self.__state.precip_nowcast_cascades = np.array( - precip_nowcast_decomp - ).swapaxes(1, 2) - self.__state.mean_nowcast_timestep = np.array( - precip_nowcast_means - ).swapaxes(1, 2) - self.__state.std_nowcast_timestep = np.array(precip_nowcast_stds).swapaxes( - 1, 2 - ) + precip_nowcast_decomp = [] + precip_nowcast_means = [] + precip_nowcast_stds = [] + + for i in range(self.__precip_nowcast.shape[0]): + precip_nowcast_decomp.append(results[i]["precip_nowcast_decomp"]) + + precip_nowcast_means.append(results[i]["precip_nowcast_means"]) + + precip_nowcast_stds.append(results[i]["precip_nowcast_stds"]) + self.__state.precip_nowcast_cascades = np.array( + precip_nowcast_decomp + ).swapaxes(1, 2) + self.__state.mean_nowcast_timestep = np.array( + precip_nowcast_means + ).swapaxes(1, 2) + self.__state.std_nowcast_timestep = np.array( + precip_nowcast_stds + ).swapaxes(1, 2) # If necessary, recompose (NWP) model forecasts self.__state.precip_models_cascades = None @@ -1790,8 +1803,9 @@ def __find_nowcast_NWP_combination(self, t): "The number of nowcast ensemble members provided is larger than the given number of ensemble members requested. n_ens_members_provided <= n_ens_members." ) - n_ens_members_max = max(n_ens_members_provided, n_model_members) + n_ens_members_max = self.__config.n_ens_members n_ens_members_min = min(n_ens_members_provided, n_model_members) + # Also make a list of the model index numbers. These indices are needed # for indexing the right climatological skill file when pysteps calculates # the blended forecast in parallel. @@ -1802,68 +1816,9 @@ def __find_nowcast_NWP_combination(self, t): else: self.__state.mapping_list_NWP_member_to_ensemble_member = [0] - # Now, repeat the nowcast ensemble members or the nwp models/members until - # it has the same amount of members as n_ens_members_max. For instance, if - # you have 10 ensemble nowcasts members and 3 NWP members, the output will - # be an ensemble of 10 members. Hence, the three NWP members are blended - # with the first three members of the nowcast (member one with member one, - # two with two, etc.), subsequently, the same NWP members are blended with - # the next three members (NWP member one with member 4, NWP member 2 with - # member 5, etc.), until 10 is reached. - if n_ens_members_min != n_ens_members_max: - if n_model_members == 1: + def repeat_precip_to_match_ensemble_size(repeats, model_type): + if model_type == "nwp": print("Repeating the NWP model for all ensemble members") - self.__state.precip_models_cascades_timestep = np.repeat( - self.__state.precip_models_cascades_timestep, - n_ens_members_max, - axis=0, - ) - self.__state.mean_models_timestep = np.repeat( - self.__state.mean_models_timestep, n_ens_members_max, axis=0 - ) - self.__state.std_models_timestep = np.repeat( - self.__state.std_models_timestep, n_ens_members_max, axis=0 - ) - self.__state.velocity_models_timestep = np.repeat( - self.__state.velocity_models_timestep, n_ens_members_max, axis=0 - ) - # For the prob. matching - self.__state.precip_models_timestep = np.repeat( - self.__state.precip_models_timestep, n_ens_members_max, axis=0 - ) - # Finally, for the model indices - self.__state.mapping_list_NWP_member_to_ensemble_member = np.repeat( - self.__state.mapping_list_NWP_member_to_ensemble_member, - n_ens_members_max, - axis=0, - ) - elif n_ens_members_provided == 1: - print("Repeating the nowcast for all ensemble members") - self.__state.precip_nowcast_cascades = np.repeat( - self.__state.precip_nowcast_cascades, - n_ens_members_max, - axis=0, - ) - self.__state.mean_nowcast_timestep = np.repeat( - self.__state.mean_nowcast_timestep, n_ens_members_max, axis=0 - ) - self.__state.std_nowcast_timestep = np.repeat( - self.__state.std_nowcast, n_ens_members_max, axis=0 - ) - self.__state.velocity_nowcast = np.repeat( - self.__state.velocity_nowcast, n_ens_members_max, axis=0 - ) - # For the prob. matching - self.__state.precip_nowcast_timestep = np.repeat( - self.__state.precip_nowcast_timestep, n_ens_members_max, axis=0 - ) - - elif n_model_members == n_ens_members_min: - print("Repeating the NWP model for all ensemble members") - repeats = [ - (n_ens_members_max + i) // n_ens_members_min - for i in range(n_ens_members_min) - ] self.__state.precip_models_cascades_timestep = np.repeat( self.__state.precip_models_cascades_timestep, repeats, @@ -1888,12 +1843,8 @@ def __find_nowcast_NWP_combination(self, t): repeats, axis=0, ) - elif n_ens_members_provided == n_ens_members_min: + if model_type == "nowcast": print("Repeating the nowcast for all ensemble members") - repeats = [ - (n_ens_members_max + i) // n_ens_members_min - for i in range(n_ens_members_min) - ] self.__state.precip_nowcast_cascades = np.repeat( self.__state.precip_nowcast_cascades, repeats, @@ -1902,17 +1853,126 @@ def __find_nowcast_NWP_combination(self, t): self.__state.mean_nowcast_timestep = np.repeat( self.__state.mean_nowcast_timestep, repeats, axis=0 ) - self.__state.std_nowcast = np.repeat( - self.__state.std_nowcast, repeats, axis=0 - ) - self.__state.velocity_nowcast = np.repeat( - self.__state.velocity_nowcast, repeats, axis=0 + self.__state.std_nowcast_timestep = np.repeat( + self.__state.std_nowcast_timestep, repeats, axis=0 ) # For the prob. matching self.__state.precip_nowcast_timestep = np.repeat( self.__state.precip_nowcast_timestep, repeats, axis=0 ) + # Now, repeat the nowcast ensemble members or the nwp models/members until + # it has the same amount of members as n_ens_members_max. For instance, if + # you have 10 ensemble nowcasts members and 3 NWP members, the output will + # be an ensemble of 10 members. Hence, the three NWP members are blended + # with the first three members of the nowcast (member one with member one, + # two with two, etc.), subsequently, the same NWP members are blended with + # the next three members (NWP member one with member 4, NWP member 2 with + # member 5, etc.), until 10 is reached. + if n_ens_members_min != n_ens_members_max: + if n_model_members == 1: + repeat_precip_to_match_ensemble_size(n_ens_members_max, "nwp") + # print("Repeating the NWP model for all ensemble members") + # self.__state.precip_models_cascades_timestep = np.repeat( + # self.__state.precip_models_cascades_timestep, + # n_ens_members_max, + # axis=0, + # ) + # self.__state.mean_models_timestep = np.repeat( + # self.__state.mean_models_timestep, n_ens_members_max, axis=0 + # ) + # self.__state.std_models_timestep = np.repeat( + # self.__state.std_models_timestep, n_ens_members_max, axis=0 + # ) + # self.__state.velocity_models_timestep = np.repeat( + # self.__state.velocity_models_timestep, n_ens_members_max, axis=0 + # ) + # # For the prob. matching + # self.__state.precip_models_timestep = np.repeat( + # self.__state.precip_models_timestep, n_ens_members_max, axis=0 + # ) + # # Finally, for the model indices + # self.__state.mapping_list_NWP_member_to_ensemble_member = np.repeat( + # self.__state.mapping_list_NWP_member_to_ensemble_member, + # n_ens_members_max, + # axis=0, + # ) + if n_ens_members_provided == 1: + repeat_precip_to_match_ensemble_size(n_ens_members_max, "nowcast") + # print("Repeating the nowcast for all ensemble members") + # self.__state.precip_nowcast_cascades = np.repeat( + # self.__state.precip_nowcast_cascades, + # n_ens_members_max, + # axis=0, + # ) + # self.__state.mean_nowcast_timestep = np.repeat( + # self.__state.mean_nowcast_timestep, n_ens_members_max, axis=0 + # ) + # self.__state.std_nowcast_timestep = np.repeat( + # self.__state.std_nowcast_timestep, n_ens_members_max, axis=0 + # ) + # # For the prob. matching + # self.__state.precip_nowcast_timestep = np.repeat( + # self.__state.precip_nowcast_timestep, n_ens_members_max, axis=0 + # ) + + if n_model_members == n_ens_members_min and n_model_members != 1: + print("Repeating the NWP model for all ensemble members") + repeats = [ + (n_ens_members_max + i) // n_ens_members_min + for i in range(n_ens_members_min) + ] + repeat_precip_to_match_ensemble_size(repeats, "nwp") + # self.__state.precip_models_cascades_timestep = np.repeat( + # self.__state.precip_models_cascades_timestep, + # repeats, + # axis=0, + # ) + # self.__state.mean_models_timestep = np.repeat( + # self.__state.mean_models_timestep, repeats, axis=0 + # ) + # self.__state.std_models_timestep = np.repeat( + # self.__state.std_models_timestep, repeats, axis=0 + # ) + # self.__state.velocity_models_timestep = np.repeat( + # self.__state.velocity_models_timestep, repeats, axis=0 + # ) + # # For the prob. matching + # self.__state.precip_models_timestep = np.repeat( + # self.__state.precip_models_timestep, repeats, axis=0 + # ) + # # Finally, for the model indices + # self.__state.mapping_list_NWP_member_to_ensemble_member = np.repeat( + # self.__state.mapping_list_NWP_member_to_ensemble_member, + # repeats, + # axis=0, + # ) + if ( + n_ens_members_provided == n_ens_members_min + and n_ens_members_provided != 1 + ): + repeat_precip_to_match_ensemble_size(repeats, "nowcast") + # print("Repeating the nowcast for all ensemble members") + # repeats = [ + # (n_ens_members_max + i) // n_ens_members_min + # for i in range(n_ens_members_min) + # ] + # self.__state.precip_nowcast_cascades = np.repeat( + # self.__state.precip_nowcast_cascades, + # repeats, + # axis=0, + # ) + # self.__state.mean_nowcast_timestep = np.repeat( + # self.__state.mean_nowcast_timestep, repeats, axis=0 + # ) + # self.__state.std_nowcast_timestep = np.repeat( + # self.__state.std_nowcast_timestep, repeats, axis=0 + # ) + # # For the prob. matching + # self.__state.precip_nowcast_timestep = np.repeat( + # self.__state.precip_nowcast_timestep, repeats, axis=0 + # ) + else: # Start with determining the maximum and mimimum number of members/models # in both input products From a5e8ea32eee7d17d98989e46e1981763cf80ec0e Mon Sep 17 00:00:00 2001 From: Joep1999 Date: Tue, 23 Sep 2025 14:01:34 +0200 Subject: [PATCH 12/18] concorporated changes ruben and small bugfix --- pysteps/blending/steps.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 9a218f935..b831e56f2 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -1850,6 +1850,11 @@ def repeat_precip_to_match_ensemble_size(repeats, model_type): repeats, axis=0, ) + self.__precip_nowcast = np.repeat( + self.__precip_nowcast, + repeats, + axis=0, + ) self.__state.mean_nowcast_timestep = np.repeat( self.__state.mean_nowcast_timestep, repeats, axis=0 ) From f690fc73e4a1fc45b8f7dd9c1ea808aa3e82eb68 Mon Sep 17 00:00:00 2001 From: Ruben Imhoff Date: Tue, 23 Sep 2025 14:23:02 +0200 Subject: [PATCH 13/18] fix: minor fixes to decomposition of external nowcast --- pysteps/blending/steps.py | 22 +++++++--------------- 1 file changed, 7 insertions(+), 15 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index f725d729b..1725da9a5 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -1129,32 +1129,24 @@ def transform_to_lagrangian(precip, i): # Decompose precomputed nowcasts and rearange them again into the required components if self.__precip_nowcast is not None: if self.__precip_nowcast.shape[0] == 1: - results = self.__decompose_member(self.__precip_nowcast[0]) + results = [self.__decompose_member(self.__precip_nowcast[0])] else: with ThreadPool(self.__config.num_workers) as pool: results = pool.map( partial(self.__decompose_member), list(self.__precip_nowcast), ) - precip_nowcast_decomp = [] - precip_nowcast_means = [] - precip_nowcast_stds = [] - for i in range(self.__precip_nowcast.shape[0]): - precip_nowcast_decomp.append(results[i]["precip_nowcast_decomp"]) - - precip_nowcast_means.append(results[i]["precip_nowcast_means"]) - - precip_nowcast_stds.append(results[i]["precip_nowcast_stds"]) self.__state.precip_nowcast_cascades = np.array( - precip_nowcast_decomp + [result["precip_nowcast_decomp"] for result in results] ).swapaxes(1, 2) self.__state.mean_nowcast_timestep = np.array( - precip_nowcast_means + [result["precip_nowcast_means"] for result in results] ).swapaxes(1, 2) - self.__state.std_nowcast_timestep = np.array(precip_nowcast_stds).swapaxes( - 1, 2 - ) + self.__state.std_nowcast_timestep = np.array( + [result["precip_nowcast_stds"] for result in results] + ).swapaxes(1, 2) + # If necessary, recompose (NWP) model forecasts self.__state.precip_models_cascades = None From 3174d7ace3fae957bb0fd75acca666104c8d8aef Mon Sep 17 00:00:00 2001 From: Ruben Imhoff Date: Tue, 23 Sep 2025 15:49:51 +0200 Subject: [PATCH 14/18] fix: make sure slicing of means and stds goes well --- pysteps/blending/steps.py | 43 +++++++++++++++++++++------------------ 1 file changed, 23 insertions(+), 20 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 3b47ec87d..e31b110a4 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -361,6 +361,8 @@ class StepsBlendingState: # Extrapolation states mean_extrapolation: np.ndarray | None = None std_extrapolation: np.ndarray | None = None + mean_nowcast: np.ndarray | None = None + std_nowcast: np.ndarray | None = None mean_nowcast_timestep: np.ndarray | None = None std_nowcast_timestep: np.ndarray | None = None rho_extrap_cascade_prev: np.ndarray | None = None @@ -635,6 +637,8 @@ def __blended_nowcast_main_loop(self): final_blended_forecast_all_members_one_timestep = [ None for _ in range(self.__config.n_ens_members) ] + self.__state.mean_nowcast_timestep = self.__state.mean_nowcast[:, :, t] + self.__state.std_nowcast_timestep = self.__state.std_nowcast[:, :, t] def worker(j): worker_state = copy(self.__state) @@ -1140,10 +1144,10 @@ def transform_to_lagrangian(precip, i): self.__state.precip_nowcast_cascades = np.array( [result["precip_nowcast_decomp"] for result in results] ).swapaxes(1, 2) - self.__state.mean_nowcast_timestep = np.array( + self.__state.mean_nowcast = np.array( [result["precip_nowcast_means"] for result in results] ).swapaxes(1, 2) - self.__state.std_nowcast_timestep = np.array( + self.__state.std_nowcast = np.array( [result["precip_nowcast_stds"] for result in results] ).swapaxes(1, 2) @@ -1834,11 +1838,11 @@ def repeat_precip_to_match_ensemble_size(repeats, model_type): repeats, axis=0, ) - self.__state.mean_nowcast_timestep = np.repeat( - self.__state.mean_nowcast_timestep, repeats, axis=0 + self.__state.mean_nowcast = np.repeat( + self.__state.mean_nowcast, repeats, axis=0 ) - self.__state.std_nowcast_timestep = np.repeat( - self.__state.std_nowcast_timestep, repeats, axis=0 + self.__state.std_nowcast = np.repeat( + self.__state.std_nowcast, repeats, axis=0 ) # For the prob. matching self.__state.precip_nowcast_timestep = np.repeat( @@ -1889,11 +1893,11 @@ def repeat_precip_to_match_ensemble_size(repeats, model_type): # n_ens_members_max, # axis=0, # ) - # self.__state.mean_nowcast_timestep = np.repeat( - # self.__state.mean_nowcast_timestep, n_ens_members_max, axis=0 + # self.__state.mean_nowcast = np.repeat( + # self.__state.mean_nowcast, n_ens_members_max, axis=0 # ) - # self.__state.std_nowcast_timestep = np.repeat( - # self.__state.std_nowcast_timestep, n_ens_members_max, axis=0 + # self.__state.std_nowcast = np.repeat( + # self.__state.std_nowcast, n_ens_members_max, axis=0 # ) # # For the prob. matching # self.__state.precip_nowcast_timestep = np.repeat( @@ -1946,11 +1950,11 @@ def repeat_precip_to_match_ensemble_size(repeats, model_type): # repeats, # axis=0, # ) - # self.__state.mean_nowcast_timestep = np.repeat( - # self.__state.mean_nowcast_timestep, repeats, axis=0 + # self.__state.mean_nowcast = np.repeat( + # self.__state.mean_nowcast, repeats, axis=0 # ) - # self.__state.std_nowcast_timestep = np.repeat( - # self.__state.std_nowcast_timestep, repeats, axis=0 + # self.__state.std_nowcast = np.repeat( + # self.__state.std_nowcast, repeats, axis=0 # ) # # For the prob. matching # self.__state.precip_nowcast_timestep = np.repeat( @@ -2564,6 +2568,7 @@ def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( worker_state.precip_extrapolated_probability_matching = np.stack( worker_state.precip_extrapolated_probability_matching ) + if len(worker_state.noise_extrapolated_decomp) > 0: if self.__config.noise_method is not None: worker_state.noise_extrapolated_decomp = np.stack( worker_state.noise_extrapolated_decomp @@ -2741,14 +2746,14 @@ def __blend_cascades(self, t_sub, j, worker_state): ) # [(extr_field, n_model_fields), n_cascade_levels, ...] means_stacked = np.concatenate( ( - worker_state.mean_nowcast_timestep[None, j, :, t_sub], + worker_state.mean_nowcast_timestep[None, j, :], worker_state.mean_models_timestep, ), axis=0, ) sigmas_stacked = np.concatenate( ( - worker_state.std_nowcast_timestep[None, j, :, t_sub], + worker_state.std_nowcast_timestep[None, j, :], worker_state.std_models_timestep, ), axis=0, @@ -2811,16 +2816,14 @@ def __blend_cascades(self, t_sub, j, worker_state): ) # [(extr_field, n_model_fields), n_cascade_levels, ...] means_stacked = np.concatenate( ( - # worker_state.mean_nowcast_timestep, - worker_state.mean_nowcast_timestep[None, j, :, t_sub], + worker_state.mean_nowcast_timestep[None, j, :], worker_state.mean_models_timestep[None, j], ), axis=0, ) sigmas_stacked = np.concatenate( ( - # worker_state.std_nowcast_timestep, - worker_state.std_nowcast_timestep[None, j, :, t_sub], + worker_state.std_nowcast_timestep[None, j, :], worker_state.std_models_timestep[None, j], ), axis=0, From a35461b7d84cc058ccd5da5aeca715d9f1b58742 Mon Sep 17 00:00:00 2001 From: Ruben Imhoff Date: Fri, 26 Sep 2025 15:48:29 +0200 Subject: [PATCH 15/18] feat: add ensemble forecast gallery example --- examples/steps_blended_forecast.py | 136 +++++++++++++++++++++++++---- 1 file changed, 120 insertions(+), 16 deletions(-) diff --git a/examples/steps_blended_forecast.py b/examples/steps_blended_forecast.py index 741679d3e..b67f6ef39 100644 --- a/examples/steps_blended_forecast.py +++ b/examples/steps_blended_forecast.py @@ -195,7 +195,7 @@ 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) + ax1 = plt.subplot(n_leadtimes, 2, n * 2 + 1) plot_precip_field( precip_forecast[0, int(leadtime / timestep) - 1, :, :], geodata=radar_metadata, @@ -204,10 +204,11 @@ colorscale="STEPS-NL", colorbar=False, ) + ax1.axis("off") # Raw NWP forecast plt.subplot(n_leadtimes, 2, n * 2 + 2) - plot_precip_field( + ax2 = plot_precip_field( nwp_precip_mmh[0, int(leadtime / timestep) - 1, :, :], geodata=nwp_metadata, title=f"NWP +{leadtime} min", @@ -215,26 +216,27 @@ colorscale="STEPS-NL", colorbar=False, ) + ax2.axis("off") plt.tight_layout() plt.show() ############################################################################### -# It is also possible to blend a deterministic external nowcast (e.g. a pre- -# made nowcast or a deterministic AI-based nowcast) with NWP using the STEPS -# algorithm. In that case, we add a `precip_nowcast` to `blending.steps.forecast`. -# By providing an external nowcast, the STEPS blending method will omit the -# autoregression and advection step for the extrapolation cascade and use the -# existing external nowcast instead (which will thus be decomposed into -# multiplicative cascades!). The weights determination and possible post- -# processings steps will remain the same. +# It is also possible to blend a deterministic or probabilistic external nowcast +# (e.g. a pre-made nowcast or a deterministic AI-based nowcast) with NWP using +# the STEPS algorithm. In that case, we add a `precip_nowcast` to +# `blending.steps.forecast`. By providing an external nowcast, the STEPS blending +# method will omit the autoregression and advection step for the extrapolation +# cascade and use the existing external nowcast instead (which will thus be +# decomposed into multiplicative cascades!). The weights determination and +# possible post-processings steps will remain the same. # # Start with creating an external nowcast # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # We go for a simple advection-only nowcast for the example, but this setup can -# be replaced with any external deterministic nowcast. +# be replaced with any external deterministic or probabilistic nowcast. extrapolate = nowcasts.get_method("extrapolation") radar_precip_to_advect = radar_precip.copy() radar_metadata_to_advect = radar_metadata.copy() @@ -258,8 +260,8 @@ ################################################################################ -# Blend the external nowcast with NWP -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# Blend the external nowcast with NWP - deterministic mode +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ precip_forecast = blending.steps.forecast( precip=radar_precip, @@ -311,7 +313,7 @@ idx = int(leadtime / timestep) - 1 # Blended nowcast - plt.subplot(n_leadtimes, 3, n * 3 + 1) + ax1 = plt.subplot(n_leadtimes, 3, n * 3 + 1) plot_precip_field( precip_forecast[0, idx, :, :], geodata=radar_metadata, @@ -320,9 +322,10 @@ colorscale="STEPS-NL", colorbar=False, ) + ax1.axis("off") # Raw extrapolated nowcast - plt.subplot(n_leadtimes, 3, n * 3 + 2) + ax2 = plt.subplot(n_leadtimes, 3, n * 3 + 2) plot_precip_field( fc_lagrangian_extrapolation_mmh[idx, :, :], geodata=radar_metadata, @@ -331,10 +334,11 @@ colorscale="STEPS-NL", colorbar=False, ) + ax2.axis("off") # Raw NWP forecast plt.subplot(n_leadtimes, 3, n * 3 + 3) - plot_precip_field( + ax3 = plot_precip_field( nwp_precip_mmh[0, idx, :, :], geodata=nwp_metadata, title=f"NWP +{leadtime} min", @@ -342,6 +346,106 @@ colorscale="STEPS-NL", colorbar=False, ) + ax3.axis("off") + + +################################################################################ +# Blend the external nowcast with NWP - ensemble mode +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +precip_forecast = blending.steps.forecast( + precip=radar_precip, + precip_nowcast=fc_lagrangian_extrapolation, + precip_models=nwp_precip, + velocity=velocity_radar, + velocity_models=velocity_nwp, + timesteps=18, + timestep=timestep, + issuetime=date_radar, + n_ens_members=5, + precip_thr=radar_metadata["threshold"], + kmperpixel=radar_metadata["xpixelsize"] / 1000.0, + noise_stddev_adj="auto", + vel_pert_method=None, + nowcasting_method="external_nowcast", + noise_method="nonparametric", + probmatching_method="cdf", + mask_method="incremental", + weights_method="bps", +) + +# Transform the data back into mm/h +precip_forecast, _ = converter(precip_forecast, radar_metadata) +radar_precip_mmh, _ = converter(radar_precip, radar_metadata) +fc_lagrangian_extrapolation_mmh, _ = converter( + fc_lagrangian_extrapolation, radar_metadata_to_advect +) +nwp_precipfc_lagrangian_extrapolation_mmh_mmh, _ = converter(nwp_precip, nwp_metadata) + + +################################################################################ +# Visualize the output +# ~~~~~~~~~~~~~~~~~~~~ + +fig = plt.figure(figsize=(5, 12)) + +leadtimes_min = [30, 60, 90, 120, 150, 180] +n_leadtimes = len(leadtimes_min) + +for n, leadtime in enumerate(leadtimes_min): + idx = int(leadtime / timestep) - 1 + + # Blended nowcast member 1 + ax1 = plt.subplot(n_leadtimes, 4, n * 4 + 1) + plot_precip_field( + precip_forecast[0, idx, :, :], + geodata=radar_metadata, + title="Blend Mem. 1", + axis="off", + colorscale="STEPS-NL", + colorbar=False, + ) + ax1.axis("off") + + # Blended nowcast member 5 + ax2 = plt.subplot(n_leadtimes, 4, n * 4 + 2) + plot_precip_field( + precip_forecast[4, idx, :, :], + geodata=radar_metadata, + title="Blend Mem. 5", + axis="off", + colorscale="STEPS-NL", + colorbar=False, + ) + ax2.axis("off") + + # Raw extrapolated nowcast + ax3 = plt.subplot(n_leadtimes, 4, n * 4 + 3) + plot_precip_field( + fc_lagrangian_extrapolation_mmh[0, idx, :, :], + geodata=radar_metadata, + title=f"NWC + {leadtime} min", + axis="off", + colorscale="STEPS-NL", + colorbar=False, + ) + ax3.axis("off") + + # Raw NWP forecast + ax4 = plt.subplot(n_leadtimes, 4, n * 4 + 4) + plot_precip_field( + nwp_precip_mmh[0, idx, :, :], + geodata=nwp_metadata, + title=f"NWP + {leadtime} min", + axis="off", + colorscale="STEPS-NL", + colorbar=False, + ) + ax4.axis("off") + +plt.show() + +print("Done.") ################################################################################ From c43fac9443b1342f645e3b3de4a3bac5a693d1f2 Mon Sep 17 00:00:00 2001 From: Joep1999 Date: Thu, 9 Oct 2025 12:26:01 +0200 Subject: [PATCH 16/18] added first test to check working of creating ensembles with deterministic data --- pysteps/tests/test_blending_steps.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/pysteps/tests/test_blending_steps.py b/pysteps/tests/test_blending_steps.py index 8ab996ce4..4e6f9aeef 100644 --- a/pysteps/tests/test_blending_steps.py +++ b/pysteps/tests/test_blending_steps.py @@ -55,17 +55,18 @@ (5, 3, 5, 6,'steps', "incremental", "cdf", False, "bps", False, 5, False, False, 80, False, None, 40), (5, 3, 5, 6,'steps', "incremental", "cdf", False, "bps", False, 5, False, False, 80, False, None, 60), #Test the externally provided nowcast - (1, 10, 1, 8,'external_nowcast', None, None, False, "spn", True, 1, False, False, 0, False, None, None), - (1, 10, 1, 8,'external_nowcast', "incremental", None, False, "bps", True, 1, False, False, 0, False, None, None), - (1, 10, 1, 8,'external_nowcast', "incremental", None, False, "spn", True, 1, False, False, 80, False, None, None), - (1, 10, 1, 8,'external_nowcast', "incremental", None, False, "bps", True, 1, True, False, 0, False, None, None), - (1, 10, 1, 8,'external_nowcast', "incremental", None, False, "spn", True, 1, False, True, 0, False, None, None), - (1, 10, 1, 8,'external_nowcast', "incremental", None, False, "bps", True, 1, True, True, 0, False, None, None), - (1, 10, 1, 8,'external_nowcast', "incremental", "cdf", False, "spn", True, 1, False, False, 0, True, None, None), - (1, 10, 1, 8,'external_nowcast', "incremental", "obs", False, "bps", True, 1, False, False, 0, False, None, None), + (1, 10, 1, 8,'external_nowcast_ens', None, None, False, "spn", True, 1, False, False, 0, False, None, None), + (1, 10, 1, 8,'external_nowcast_ens', "incremental", None, False, "bps", True, 1, False, False, 0, False, None, None), + (1, 10, 1, 8,'external_nowcast_ens', "incremental", None, False, "spn", True, 1, False, False, 80, False, None, None), + (1, 10, 1, 8,'external_nowcast_ens', "incremental", None, False, "bps", True, 1, True, False, 0, False, None, None), + (1, 10, 1, 8,'external_nowcast_ens', "incremental", None, False, "spn", True, 1, False, True, 0, False, None, None), + (1, 10, 1, 8,'external_nowcast_ens', "incremental", None, False, "bps", True, 1, True, True, 0, False, None, None), + (1, 10, 1, 8,'external_nowcast_ens', "incremental", "cdf", False, "spn", True, 1, False, False, 0, True, None, None), + (1, 10, 1, 8,'external_nowcast_ens', "incremental", "obs", False, "bps", True, 1, False, False, 0, False, None, None), (5, 10, 5, 8,'external_nowcast_det', "incremental", None, False, "spn", True, 5, False, False, 0, False, None, None), (5, 10, 5, 8,'external_nowcast_ens', "incremental", None, False, "spn", True, 5, False, False, 0, False, None, None), - (1, 10, 5, 8,'external_nowcast_det', "incremental", None, False, "spn", True, 5, False, False, 0, False, None, None), + (1, 10, 5, 8,'external_nowcast_ens', "incremental", None, False, "spn", True, 5, False, False, 0, False, None, None), + (1, 10, 1, 8,'external_nowcast_det', "incremental", "cdf", False, "spn", True, 5, False, False, 0, False, None, None) ] # fmt:on @@ -214,7 +215,7 @@ def test_steps_blending( n_ens_member, i, 30:165, 41 + 1 * (i + 1) * n_ens_member ] = 0.1 - elif nowcasting_method[:16] == "external_nowcast": + elif nowcasting_method == "external_nowcast_det": nowcasting_method = "external_nowcast" for i in range(precip_nowcast.shape[1]): precip_nowcast[0, i, 30:165, 30 + 1 * i] = 0.1 From 7e04001b1e8769ebbf4f35972b3387ffdabb7f73 Mon Sep 17 00:00:00 2001 From: Joep1999 Date: Thu, 9 Oct 2025 15:21:36 +0200 Subject: [PATCH 17/18] small fixes tests and blending --- pysteps/blending/steps.py | 5 +++-- pysteps/tests/test_blending_steps.py | 20 +++++++++++++------- 2 files changed, 16 insertions(+), 9 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index e31b110a4..7cde87fd4 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -637,8 +637,9 @@ def __blended_nowcast_main_loop(self): final_blended_forecast_all_members_one_timestep = [ None for _ in range(self.__config.n_ens_members) ] - self.__state.mean_nowcast_timestep = self.__state.mean_nowcast[:, :, t] - self.__state.std_nowcast_timestep = self.__state.std_nowcast[:, :, t] + if self.__config.nowcasting_method == "external_nowcast": + self.__state.mean_nowcast_timestep = self.__state.mean_nowcast[:, :, t] + self.__state.std_nowcast_timestep = self.__state.std_nowcast[:, :, t] def worker(j): worker_state = copy(self.__state) diff --git a/pysteps/tests/test_blending_steps.py b/pysteps/tests/test_blending_steps.py index 4e6f9aeef..95fbb7e99 100644 --- a/pysteps/tests/test_blending_steps.py +++ b/pysteps/tests/test_blending_steps.py @@ -66,8 +66,10 @@ (5, 10, 5, 8,'external_nowcast_det', "incremental", None, False, "spn", True, 5, False, False, 0, False, None, None), (5, 10, 5, 8,'external_nowcast_ens', "incremental", None, False, "spn", True, 5, False, False, 0, False, None, None), (1, 10, 5, 8,'external_nowcast_ens', "incremental", None, False, "spn", True, 5, False, False, 0, False, None, None), - (1, 10, 1, 8,'external_nowcast_det', "incremental", "cdf", False, "spn", True, 5, False, False, 0, False, None, None) + (1, 10, 1, 8,'external_nowcast_det', "incremental", "cdf", False, "bps", True, 5, False, False, 0, False, None, None), + (1, 10, 1, 8,'external_nowcast_det', "incremental", "obs", False, "spn", True, 5, False, False, 0, False, None, None) ] + # fmt:on steps_arg_names = ( @@ -391,9 +393,7 @@ def test_steps_blending( ) assert precip_forecast.ndim == 4, "Wrong amount of dimensions in forecast output" - assert ( - precip_forecast.shape[0] == expected_n_ens_members - ), "Wrong amount of output ensemble members in forecast output" + assert ( precip_forecast.shape[1] == n_timesteps ), "Wrong amount of output time steps in forecast output" @@ -404,9 +404,15 @@ def test_steps_blending( assert ( precip_forecast.ndim == 4 ), "Wrong amount of dimensions in converted forecast output" - assert ( - precip_forecast.shape[0] == expected_n_ens_members - ), "Wrong amount of output ensemble members in converted forecast output" + if nowcasting_method != "external_nowcast": + assert ( + precip_forecast.shape[0] == expected_n_ens_members + ), "Wrong amount of output ensemble members in forecast output" + + assert ( + precip_forecast.shape[0] == expected_n_ens_members + ), "Wrong amount of output ensemble members in converted forecast output" + assert ( precip_forecast.shape[1] == n_timesteps ), "Wrong amount of output time steps in converted forecast output" From 3b256e934aea203c4adeab8e4ac5d067c4951d8b Mon Sep 17 00:00:00 2001 From: Joep1999 Date: Thu, 9 Oct 2025 16:40:15 +0200 Subject: [PATCH 18/18] fixed a mistake in the small fixes --- pysteps/tests/test_blending_steps.py | 40 +++++++++++++++------------- 1 file changed, 21 insertions(+), 19 deletions(-) diff --git a/pysteps/tests/test_blending_steps.py b/pysteps/tests/test_blending_steps.py index 95fbb7e99..f914638f7 100644 --- a/pysteps/tests/test_blending_steps.py +++ b/pysteps/tests/test_blending_steps.py @@ -55,19 +55,19 @@ (5, 3, 5, 6,'steps', "incremental", "cdf", False, "bps", False, 5, False, False, 80, False, None, 40), (5, 3, 5, 6,'steps', "incremental", "cdf", False, "bps", False, 5, False, False, 80, False, None, 60), #Test the externally provided nowcast - (1, 10, 1, 8,'external_nowcast_ens', None, None, False, "spn", True, 1, False, False, 0, False, None, None), - (1, 10, 1, 8,'external_nowcast_ens', "incremental", None, False, "bps", True, 1, False, False, 0, False, None, None), - (1, 10, 1, 8,'external_nowcast_ens', "incremental", None, False, "spn", True, 1, False, False, 80, False, None, None), - (1, 10, 1, 8,'external_nowcast_ens', "incremental", None, False, "bps", True, 1, True, False, 0, False, None, None), - (1, 10, 1, 8,'external_nowcast_ens', "incremental", None, False, "spn", True, 1, False, True, 0, False, None, None), - (1, 10, 1, 8,'external_nowcast_ens', "incremental", None, False, "bps", True, 1, True, True, 0, False, None, None), - (1, 10, 1, 8,'external_nowcast_ens', "incremental", "cdf", False, "spn", True, 1, False, False, 0, True, None, None), - (1, 10, 1, 8,'external_nowcast_ens', "incremental", "obs", False, "bps", True, 1, False, False, 0, False, None, None), - (5, 10, 5, 8,'external_nowcast_det', "incremental", None, False, "spn", True, 5, False, False, 0, False, None, None), + (1, 10, 1, 8,'external_nowcast_det', None, None, False, "spn", True, 1, False, False, 0, False, None, None), + (1, 10, 1, 8,'external_nowcast_det', "incremental", None, False, "bps", True, 1, False, False, 0, False, None, None), + (1, 10, 1, 8,'external_nowcast_det', "incremental", None, False, "spn", True, 1, False, False, 80, False, None, None), + (1, 10, 1, 8,'external_nowcast_det', "incremental", None, False, "bps", True, 1, True, False, 0, False, None, None), + (1, 10, 1, 8,'external_nowcast_det', "incremental", None, False, "spn", True, 1, False, True, 0, False, None, None), + (1, 10, 1, 8,'external_nowcast_det', "incremental", None, False, "bps", True, 1, True, True, 0, False, None, None), + (1, 10, 1, 8,'external_nowcast_det', "incremental", "cdf", False, "spn", True, 1, False, False, 0, True, None, None), + (1, 10, 1, 8,'external_nowcast_det', "incremental", "obs", False, "bps", True, 1, False, False, 0, False, None, None), + (5, 10, 5, 8,'external_nowcast_ens', "incremental", None, False, "spn", True, 5, False, False, 0, False, None, None), (5, 10, 5, 8,'external_nowcast_ens', "incremental", None, False, "spn", True, 5, False, False, 0, False, None, None), (1, 10, 5, 8,'external_nowcast_ens', "incremental", None, False, "spn", True, 5, False, False, 0, False, None, None), - (1, 10, 1, 8,'external_nowcast_det', "incremental", "cdf", False, "bps", True, 5, False, False, 0, False, None, None), - (1, 10, 1, 8,'external_nowcast_det', "incremental", "obs", False, "spn", True, 5, False, False, 0, False, None, None) + (1, 10, 1, 8,'external_nowcast_ens', "incremental", "cdf", False, "bps", True, 5, False, False, 0, False, None, None), + (5, 10, 1, 8,'external_nowcast_ens', "incremental", "obs", False, "spn", True, 5, False, False, 0, False, None, None) ] # fmt:on @@ -216,6 +216,8 @@ def test_steps_blending( precip_nowcast[ n_ens_member, i, 30:165, 41 + 1 * (i + 1) * n_ens_member ] = 0.1 + if n_ens_members < expected_n_ens_members: + n_ens_members = expected_n_ens_members elif nowcasting_method == "external_nowcast_det": nowcasting_method = "external_nowcast" @@ -394,6 +396,10 @@ def test_steps_blending( assert precip_forecast.ndim == 4, "Wrong amount of dimensions in forecast output" + assert ( + precip_forecast.shape[0] == expected_n_ens_members + ), "Wrong amount of output ensemble members in forecast output" + assert ( precip_forecast.shape[1] == n_timesteps ), "Wrong amount of output time steps in forecast output" @@ -404,14 +410,10 @@ def test_steps_blending( assert ( precip_forecast.ndim == 4 ), "Wrong amount of dimensions in converted forecast output" - if nowcasting_method != "external_nowcast": - assert ( - precip_forecast.shape[0] == expected_n_ens_members - ), "Wrong amount of output ensemble members in forecast output" - - assert ( - precip_forecast.shape[0] == expected_n_ens_members - ), "Wrong amount of output ensemble members in converted forecast output" + + assert ( + precip_forecast.shape[0] == expected_n_ens_members + ), "Wrong amount of output ensemble members in converted forecast output" assert ( precip_forecast.shape[1] == n_timesteps