diff --git a/pysteps/blending/clim.py b/pysteps/blending/clim.py index 82656b27..886fcb11 100644 --- a/pysteps/blending/clim.py +++ b/pysteps/blending/clim.py @@ -193,8 +193,13 @@ def calc_clim_skill( past_skill = np.array(None) # check if there is enough data to compute the climatological skill if not past_skill.any(): + print("WARNING: Past skill file is empty, using default BPS2006 skill") return get_default_skill(n_cascade_levels, n_models) elif past_skill.shape[0] < window_length: + print( + f"WARNING: Past skill file has fewer days ({past_skill.shape[0]})" + f"than expected ({window_length}). Using default BPS2006 skill" + ) return get_default_skill(n_cascade_levels, n_models) # reduce window if necessary else: diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index fc1530f0..9080eb4a 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -416,8 +416,6 @@ class StepsBlendingState: # Final outputs final_blended_forecast: np.ndarray | None = None final_blended_forecast_non_perturbed: np.ndarray | None = None - weights: np.ndarray | None = None - weights_model_only: np.ndarray | None = None # Timing and indexing time_prev_timestep: list[float] | None = None @@ -426,10 +424,24 @@ class StepsBlendingState: is_nowcast_time_step: bool | None = None subtimestep_index: int | None = None - # Weights used for blending + # Weights used for blending. `weights` / `weights_model_only` are per-worker + # scratch (the current member's weights, used by the downstream blends). + # `weights_per_member` / `weights_model_only_per_member` persist each member's + # weights across timesteps so members do not overwrite one another and the + # full-NWP smoothing reads each member's own previous-timestep weights. + + # Current working copy: full weights (extrapolation + NWP + noise) weights: np.ndarray | None = None + + # Current working copy: Extrapolation + radar removed (for when there is no radar) weights_model_only: np.ndarray | None = None + # List index by j that survives over timesteps (full weights) + weights_per_member: list[Any] | None = None + + # List index by j that survives over timesteps (Extrapolation + radar removed) + weights_model_only_per_member: list[Any] | None = None + # This is stores here as well because this is changed during the forecast loop and thus no longer part of the config extrapolation_kwargs: dict[str, Any] = field(default_factory=dict) @@ -665,7 +677,7 @@ def __blended_nowcast_main_loop(self): def worker(j): worker_state = copy(self.__state) self.__determine_NWP_skill_for_next_timestep(t, j, worker_state) - self.__determine_weights_per_component(t, worker_state) + self.__determine_weights_per_component(t, j, worker_state) self.__regress_extrapolation_and_noise_cascades(j, worker_state, t) self.__perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep( t, j, worker_state @@ -1629,6 +1641,14 @@ def __prepare_forecast_loop(self): self.__state.previous_displacement_prob_matching = np.stack( [None for j in range(self.__config.n_ens_members)] ) + # Per-member blending weights, persisted across timesteps so members do + # not overwrite one another (see StepsBlendingState). + self.__state.weights_per_member = [ + None for j in range(self.__config.n_ens_members) + ] + self.__state.weights_model_only_per_member = [ + None for j in range(self.__config.n_ens_members) + ] self.__state.final_blended_forecast = [ [] for j in range(self.__config.n_ens_members) ] @@ -2128,11 +2148,16 @@ def __determine_NWP_skill_for_next_timestep(self, t, j, worker_state): axis=0, ) - def __determine_weights_per_component(self, t, worker_state): + def __determine_weights_per_component(self, t, j, worker_state): """ Compute blending weights for each component based on the selected method ('bps' or 'spn'). Weights are determined for both full blending and model-only scenarios, accounting for correlations and covariance. + + The resulting weights are stored per ensemble member ``j`` in + ``weights_per_member`` / ``weights_model_only_per_member`` so members do + not overwrite one another and the full-NWP smoothing reads each member's + own previous-timestep weights. """ start_smoothing_to_final_weights = False if self.__config.timestep_start_full_nwp_weight is not None: @@ -2149,7 +2174,7 @@ def __determine_weights_per_component(self, t, worker_state): ) else: worker_state.weights = calculate_end_weights( - previous_weights=self.__state.weights, + previous_weights=worker_state.weights_per_member[j], timestep=t, n_timesteps=self.__timesteps[-1], start_full_nwp_weight=self.__config.timestep_start_full_nwp_weight, @@ -2212,7 +2237,7 @@ def __determine_weights_per_component(self, t, worker_state): ) elif start_smoothing_to_final_weights: worker_state.weights_model_only = calculate_end_weights( - previous_weights=self.__state.weights_model_only, + previous_weights=worker_state.weights_model_only_per_member[j], timestep=t, n_timesteps=self.__timesteps[-1], start_full_nwp_weight=self.__config.timestep_start_full_nwp_weight, @@ -2223,8 +2248,11 @@ def __determine_weights_per_component(self, t, worker_state): "Unknown weights method %s: must be 'bps' or 'spn'" % self.__config.weights_method ) - self.__state.weights = worker_state.weights - self.__state.weights_model_only = worker_state.weights_model_only + # Persist this member's weights without overwriting other members. The + # shallow copy shares the same list object, so index assignment touches a + # disjoint element per member (mirrors previous_displacement[j]). + worker_state.weights_per_member[j] = worker_state.weights + worker_state.weights_model_only_per_member[j] = worker_state.weights_model_only def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): """ @@ -2921,7 +2949,8 @@ def __blend_cascades(self, t_sub, j, worker_state): covariance=covariance_nwp_models, ) - self.__state.weights = worker_state.weights + # Persist this member's (spn) weights without overwriting other members. + worker_state.weights_per_member[j] = worker_state.weights # 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.