From cbd0d06679e41de127ea1100c57050d521e27123 Mon Sep 17 00:00:00 2001 From: FelixE91 Date: Thu, 26 Mar 2026 18:01:25 +0100 Subject: [PATCH 1/6] Fix precip-noprecip input array blending failure --- pysteps/blending/steps.py | 63 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 62 insertions(+), 1 deletion(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 161c4454..9011408d 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -1196,12 +1196,21 @@ def transform_to_lagrangian(precip, i): self.__precip_models = np.stack(temp_precip_models) # Check for zero input fields in the radar, nowcast and NWP data. + # decide based on latest input file only self.__params.zero_precip_radar = check_norain( - self.__precip, + self.__precip[-1], self.__params.precip_threshold, self.__config.norain_threshold, self.__params.noise_kwargs["win_fun"], ) + # Set all inputs to 0 (also previous times) + # not strictly necessary but to be consistent with checking only + # latest observation rain field to set zero_precip_radar + if self.__params.zero_precip_radar: + self.__precip = np.where(np.isfinite(self.__precip), + np.ones(self.__precip.shape)*np.nanmin(self.__precip), + self.__precip + ) # The norain fraction threshold used for nwp is the default value of 0.0, # since nwp does not suffer from clutter. @@ -1481,6 +1490,7 @@ def __estimate_ar_parameters_radar(self): # Finally, transpose GAMMA to ensure that the shape is the same as np.empty((n_cascade_levels, ar_order)) GAMMA = GAMMA.transpose() + if len(GAMMA.shape)==1: GAMMA = GAMMA.reshape((GAMMA.size,1)) assert GAMMA.shape == ( self.__config.n_cascade_levels, self.__config.ar_order, @@ -3625,6 +3635,11 @@ def forecast( turns out to be a warranted functionality. """ + # Check the input precip and ar_order to be consistent + # zero-precip in previous time steps has to be removed + # (constant field causes autoregression to fail) + precip, ar_order = check_previous_radar_obs(precip,ar_order) + blending_config = StepsBlendingConfig( n_ens_members=n_ens_members, n_cascade_levels=n_cascade_levels, @@ -3686,6 +3701,52 @@ def forecast( return forecast_steps_nowcast +# TODO: Where does this piece of code best fit: in utils or inside the class? +def check_previous_radar_obs(precip,ar_order): + """Check all radar time steps and remove zero precipitation time steps + + Parameters + ---------- + precip : array-like + Array of shape (ar_order+1,m,n) containing the input precipitation fields + ordered by timestamp from oldest to newest. The time steps between the + inputs are assumed to be regular. + ar_order : int + The order of the autoregressive model to use. Must be >= 1. + + Returns + ------- + precip : numpy array + Array of shape (ar_order+1,m,n) containing the modified array with + input precipitation fields ordered by timestamp from oldest to newest. + The time steps between the inputs are assumed to be regular. + ar_order : int + The order of the autoregressive model to use. Must be >= 1. + Adapted to match with precip.shape equal (ar_order+1,m,n). + """ + # Quick check radar input - at least 2 time steps + assert precip.shape[0] >= 2 + + # Check all time steps for zero-precip (constant field, minimum==maximum) + zero_precip = [np.nanmin(obs)==np.nanmax(obs) for obs in precip] + if zero_precip[-1] or ~np.any(zero_precip): + # Unchanged if no rain in latest time step -> will be processed as zero_precip_radar=True + # or Unchanged if all time steps contain rain + return precip, ar_order + elif zero_precip[-2]: + # try to use a previous time step + if not np.all(zero_precip[:-2]): + # find latest non-zero precip + # ATTENTION: This changes the time between precip[-2] and precip[-1] from initial 5min to a longer period + prev = np.arange(len(zero_precip[:-2]))[~zero_precip][-1] + return precip[[prev,-1]], 1 + raise ValueError("Precipitation in latest but no previous time step. Not possible to calculate autoregression.") + else: + # Keep the latest time steps that do all contain precip + precip = precip[np.max(np.arange(len(zero_precip))[zero_precip])+1:].copy() + return precip, min(ar_order,precip.shape[0]-1) + + # TODO: Where does this piece of code best fit: in utils or inside the class? def calculate_ratios(correlations): """Calculate explained variance ratios from correlation. From 4d01c01ddb5b37a6d0616102c73ec4c5ca11891e Mon Sep 17 00:00:00 2001 From: FelixE91 Date: Thu, 26 Mar 2026 18:02:10 +0100 Subject: [PATCH 2/6] Add tests for blending with noprecip and precip input time steps --- pysteps/tests/test_blending_steps.py | 152 +++++++++++++++++++++++++++ 1 file changed, 152 insertions(+) diff --git a/pysteps/tests/test_blending_steps.py b/pysteps/tests/test_blending_steps.py index fe9e2863..bd5d70b2 100644 --- a/pysteps/tests/test_blending_steps.py +++ b/pysteps/tests/test_blending_steps.py @@ -426,3 +426,155 @@ def test_steps_blending( assert ( precip_forecast.shape[1] == n_timesteps ), "Wrong amount of output time steps in converted forecast output" + + ### Add another test for observations that are 0 at last time step and non-zero before that + radar_precip_zero_latest = radar_precip.copy() + radar_precip_zero_latest[-1] = metadata["zerovalue"] + + ### + # The blending + ### + precip_forecast = blending.steps.forecast( + precip=radar_precip_zero_latest, + precip_models=nwp_precip_decomp, + velocity=radar_velocity, + velocity_models=nwp_velocity, + timesteps=timesteps, + timestep=5.0, + issuetime=datetime.datetime.strptime("202112012355", "%Y%m%d%H%M"), + n_ens_members=n_ens_members, + n_cascade_levels=n_cascade_levels, + blend_nwp_members=blend_nwp_members, + precip_thr=metadata["threshold"], + kmperpixel=1.0, + extrap_method="semilagrangian", + decomp_method="fft", + bandpass_filter_method="gaussian", + noise_method="nonparametric", + noise_stddev_adj="auto", + ar_order=2, + vel_pert_method=vel_pert_method, + weights_method=weights_method, + timestep_start_full_nwp_weight=timestep_start_full_nwp_weight, + conditional=False, + probmatching_method=probmatching_method, + mask_method=mask_method, + resample_distribution=resample_distribution, + smooth_radar_mask_range=smooth_radar_mask_range, + callback=None, + return_output=True, + seed=None, + num_workers=1, + fft_method="numpy", + domain="spatial", + outdir_path_skill=outdir_path_skill, + extrap_kwargs=None, + filter_kwargs=None, + noise_kwargs=None, + vel_pert_kwargs=None, + clim_kwargs=clim_kwargs, + mask_kwargs=mask_kwargs, + measure_time=False, + ) + + 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" + + # Transform the data back into mm/h + precip_forecast, _ = converter(precip_forecast, metadata) + + 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" + + assert ( + precip_forecast.shape[1] == n_timesteps + ), "Wrong amount of output time steps in converted forecast output" + + ### Add the case where the precip observation is non-zero for the latest time step but 0 before + radar_precip_zero_previous = radar_precip.copy() + radar_precip_zero_previous[:-1] = metadata["zerovalue"] + + ### + # The blending + ### + precip_forecast = blending.steps.forecast( + precip=radar_precip_zero_previous, + precip_models=nwp_precip_decomp, + velocity=radar_velocity, + velocity_models=nwp_velocity, + timesteps=timesteps, + timestep=5.0, + issuetime=datetime.datetime.strptime("202112012355", "%Y%m%d%H%M"), + n_ens_members=n_ens_members, + n_cascade_levels=n_cascade_levels, + blend_nwp_members=blend_nwp_members, + precip_thr=metadata["threshold"], + kmperpixel=1.0, + extrap_method="semilagrangian", + decomp_method="fft", + bandpass_filter_method="gaussian", + noise_method="nonparametric", + noise_stddev_adj="auto", + ar_order=2, + vel_pert_method=vel_pert_method, + weights_method=weights_method, + timestep_start_full_nwp_weight=timestep_start_full_nwp_weight, + conditional=False, + probmatching_method=probmatching_method, + mask_method=mask_method, + resample_distribution=resample_distribution, + smooth_radar_mask_range=smooth_radar_mask_range, + callback=None, + return_output=True, + seed=None, + num_workers=1, + fft_method="numpy", + domain="spatial", + outdir_path_skill=outdir_path_skill, + extrap_kwargs=None, + filter_kwargs=None, + noise_kwargs=None, + vel_pert_kwargs=None, + clim_kwargs=clim_kwargs, + mask_kwargs=mask_kwargs, + measure_time=False, + ) + + 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" + + # Transform the data back into mm/h + precip_forecast, _ = converter(precip_forecast, metadata) + + 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" + + assert ( + precip_forecast.shape[1] == n_timesteps + ), "Wrong amount of output time steps in converted forecast output" + + From 3948d5d9a1ada48251cbaa0654c76f636f6253f4 Mon Sep 17 00:00:00 2001 From: FelixE91 Date: Thu, 26 Mar 2026 18:46:07 +0100 Subject: [PATCH 3/6] fix prev in check_previous_radar_obs() --- pysteps/blending/steps.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 9011408d..baeff660 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -3738,7 +3738,7 @@ def check_previous_radar_obs(precip,ar_order): if not np.all(zero_precip[:-2]): # find latest non-zero precip # ATTENTION: This changes the time between precip[-2] and precip[-1] from initial 5min to a longer period - prev = np.arange(len(zero_precip[:-2]))[~zero_precip][-1] + prev = np.arange(len(zero_precip[:-2]))[~np.array(zero_precip[:-2])][-1] return precip[[prev,-1]], 1 raise ValueError("Precipitation in latest but no previous time step. Not possible to calculate autoregression.") else: From 517076810a48475a32a4764155a53de46083083c Mon Sep 17 00:00:00 2001 From: Lesley De Cruz Date: Fri, 27 Mar 2026 13:43:34 +0100 Subject: [PATCH 4/6] Improve tests for blended forecast Reduce overall number of tests and refactor code Add dedicated AR order parameterization --- pysteps/tests/test_blending_steps.py | 306 ++++++++++++--------------- 1 file changed, 140 insertions(+), 166 deletions(-) diff --git a/pysteps/tests/test_blending_steps.py b/pysteps/tests/test_blending_steps.py index bd5d70b2..b2712c68 100644 --- a/pysteps/tests/test_blending_steps.py +++ b/pysteps/tests/test_blending_steps.py @@ -77,6 +77,35 @@ # fmt:on + +def run_and_assert_forecast( + precip, forecast_kwargs, expected_n_ens_members, n_timesteps, converter, metadata +): + """Run a blended nowcast and assert the output has the expected shape.""" + precip_forecast = blending.steps.forecast(precip=precip, **forecast_kwargs) + + 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" + + # Transform the data back into mm/h + precip_forecast, _ = converter(precip_forecast, metadata) + + 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" + assert ( + precip_forecast.shape[1] == n_timesteps + ), "Wrong amount of output time steps in converted forecast output" + + steps_arg_names = ( "n_models", "timesteps", @@ -357,10 +386,9 @@ def test_steps_blending( assert nwp_velocity.ndim == 5, "nwp_velocity must be a five-dimensional array" ### - # The blending + # Shared forecast kwargs ### - precip_forecast = blending.steps.forecast( - precip=radar_precip, + forecast_kwargs = dict( precip_models=nwp_precip_decomp, velocity=radar_velocity, velocity_models=nwp_velocity, @@ -402,179 +430,125 @@ def test_steps_blending( measure_time=False, ) - 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" - - # Transform the data back into mm/h - precip_forecast, _ = converter(precip_forecast, metadata) - - 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" - - assert ( - precip_forecast.shape[1] == n_timesteps - ), "Wrong amount of output time steps in converted forecast output" - - ### Add another test for observations that are 0 at last time step and non-zero before that - radar_precip_zero_latest = radar_precip.copy() - radar_precip_zero_latest[-1] = metadata["zerovalue"] - ### # The blending ### - precip_forecast = blending.steps.forecast( - precip=radar_precip_zero_latest, - precip_models=nwp_precip_decomp, - velocity=radar_velocity, - velocity_models=nwp_velocity, - timesteps=timesteps, - timestep=5.0, - issuetime=datetime.datetime.strptime("202112012355", "%Y%m%d%H%M"), - n_ens_members=n_ens_members, - n_cascade_levels=n_cascade_levels, - blend_nwp_members=blend_nwp_members, - precip_thr=metadata["threshold"], - kmperpixel=1.0, - extrap_method="semilagrangian", - decomp_method="fft", - bandpass_filter_method="gaussian", - noise_method="nonparametric", - noise_stddev_adj="auto", - ar_order=2, - vel_pert_method=vel_pert_method, - weights_method=weights_method, - timestep_start_full_nwp_weight=timestep_start_full_nwp_weight, - conditional=False, - probmatching_method=probmatching_method, - mask_method=mask_method, - resample_distribution=resample_distribution, - smooth_radar_mask_range=smooth_radar_mask_range, - callback=None, - return_output=True, - seed=None, - num_workers=1, - fft_method="numpy", - domain="spatial", - outdir_path_skill=outdir_path_skill, - extrap_kwargs=None, - filter_kwargs=None, - noise_kwargs=None, - vel_pert_kwargs=None, - clim_kwargs=clim_kwargs, - mask_kwargs=mask_kwargs, - measure_time=False, + # Test with full radar data + run_and_assert_forecast( + radar_precip, + forecast_kwargs, + expected_n_ens_members, + n_timesteps, + converter, + metadata, ) - 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" - - # Transform the data back into mm/h - precip_forecast, _ = converter(precip_forecast, metadata) - - 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" - - assert ( - precip_forecast.shape[1] == n_timesteps - ), "Wrong amount of output time steps in converted forecast output" - - ### Add the case where the precip observation is non-zero for the latest time step but 0 before - radar_precip_zero_previous = radar_precip.copy() - radar_precip_zero_previous[:-1] = metadata["zerovalue"] +@pytest.mark.parametrize("ar_order", [1, 2]) +def test_steps_blending_partial_zero_radar(ar_order): + """Test that a forecast succeeds when only the latest radar frame has + precipitation (initiating cell corner case that produces NaN autocorrelations + for the earlier, all-zero cascade levels).""" + pytest.importorskip("cv2") - ### - # The blending - ### - precip_forecast = blending.steps.forecast( - precip=radar_precip_zero_previous, - precip_models=nwp_precip_decomp, - velocity=radar_velocity, - velocity_models=nwp_velocity, - timesteps=timesteps, - timestep=5.0, - issuetime=datetime.datetime.strptime("202112012355", "%Y%m%d%H%M"), - n_ens_members=n_ens_members, - n_cascade_levels=n_cascade_levels, - blend_nwp_members=blend_nwp_members, - precip_thr=metadata["threshold"], - kmperpixel=1.0, - extrap_method="semilagrangian", - decomp_method="fft", - bandpass_filter_method="gaussian", - noise_method="nonparametric", - noise_stddev_adj="auto", - ar_order=2, - vel_pert_method=vel_pert_method, - weights_method=weights_method, - timestep_start_full_nwp_weight=timestep_start_full_nwp_weight, - conditional=False, - probmatching_method=probmatching_method, - mask_method=mask_method, - resample_distribution=resample_distribution, - smooth_radar_mask_range=smooth_radar_mask_range, - callback=None, - return_output=True, - seed=None, - num_workers=1, - fft_method="numpy", - domain="spatial", - outdir_path_skill=outdir_path_skill, - extrap_kwargs=None, - filter_kwargs=None, - noise_kwargs=None, - vel_pert_kwargs=None, - clim_kwargs=clim_kwargs, - mask_kwargs=mask_kwargs, - measure_time=False, + n_timesteps = 3 + metadata = dict( + unit="mm", + transformation="dB", + accutime=5.0, + transform="dB", + zerovalue=0.0, + threshold=0.01, + zr_a=200.0, + zr_b=1.6, ) - 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" - - # Transform the data back into mm/h - precip_forecast, _ = converter(precip_forecast, metadata) - - assert ( - precip_forecast.ndim == 4 - ), "Wrong amount of dimensions in converted forecast output" + # Build minimal NWP data (1 model, 4 time steps) + nwp_precip = np.zeros((1, n_timesteps + 1, 200, 200)) + for i in range(nwp_precip.shape[1]): + nwp_precip[0, i, 30:185, 32 + i] = 1.0 + nwp_precip[0, i, 30:185, 33 + i] = 5.0 + nwp_precip[0, i, 30:185, 34 + i] = 5.0 + nwp_precip[0, i, 30:185, 35 + i] = 4.5 - assert ( - precip_forecast.shape[0] == expected_n_ens_members - ), "Wrong amount of output ensemble members in converted forecast output" + # Build radar data: only the latest (most recent) frame has precipitation + radar_precip = np.zeros((3, 200, 200)) + radar_precip[2, 30:155, 32] = 1.0 + radar_precip[2, 30:155, 33] = 5.0 + radar_precip[2, 30:155, 34] = 5.0 + radar_precip[2, 30:155, 35] = 4.5 - assert ( - precip_forecast.shape[1] == n_timesteps - ), "Wrong amount of output time steps in converted forecast output" + # Threshold, convert and transform + radar_precip[radar_precip < metadata["threshold"]] = 0.0 + nwp_precip[nwp_precip < metadata["threshold"]] = 0.0 + converter = pysteps.utils.get_method("mm/h") + radar_precip, _ = converter(radar_precip, metadata) + nwp_precip, metadata = converter(nwp_precip, metadata) + transformer = pysteps.utils.get_method(metadata["transformation"]) + radar_precip, _ = transformer(radar_precip, metadata) + nwp_precip, metadata = transformer(nwp_precip, metadata) + radar_precip[~np.isfinite(radar_precip)] = metadata["zerovalue"] + nwp_precip[~np.isfinite(nwp_precip)] = metadata["zerovalue"] + # Decompose NWP + n_cascade_levels = 6 + decomp_method, _ = cascade.get_method("fft") + bp_filter = cascade.get_method("gaussian")(radar_precip.shape[1:], n_cascade_levels) + nwp_precip_decomp = nwp_precip.copy() + # Velocity fields + oflow_method = pysteps.motion.get_method("lucaskanade") + radar_velocity = oflow_method(radar_precip) + _V_NWP = [ + oflow_method(nwp_precip[0, t - 1 : t + 1, :]) + for t in range(1, nwp_precip.shape[1]) + ] + nwp_velocity = np.stack([np.insert(_V_NWP, 0, _V_NWP[0], axis=0)]) + + run_and_assert_forecast( + radar_precip, + dict( + precip_models=nwp_precip_decomp, + velocity=radar_velocity, + velocity_models=nwp_velocity, + timesteps=n_timesteps, + timestep=5.0, + issuetime=datetime.datetime.strptime("202112012355", "%Y%m%d%H%M"), + n_ens_members=2, + n_cascade_levels=n_cascade_levels, + blend_nwp_members=False, + precip_thr=metadata["threshold"], + kmperpixel=1.0, + extrap_method="semilagrangian", + decomp_method="fft", + bandpass_filter_method="gaussian", + noise_method="nonparametric", + noise_stddev_adj="auto", + ar_order=ar_order, + vel_pert_method=None, + weights_method="bps", + conditional=False, + probmatching_method=None, + mask_method="incremental", + resample_distribution=False, + smooth_radar_mask_range=0, + callback=None, + return_output=True, + seed=42, + num_workers=1, + fft_method="numpy", + domain="spatial", + outdir_path_skill="./tmp/", + extrap_kwargs=None, + filter_kwargs=None, + noise_kwargs=None, + vel_pert_kwargs=None, + clim_kwargs=None, + mask_kwargs=None, + measure_time=False, + ), + expected_n_ens_members=2, + n_timesteps=n_timesteps, + converter=converter, + metadata=metadata, + ) From f3355c1b9224c78143e84e9d72f0e9d77267ae0a Mon Sep 17 00:00:00 2001 From: Lesley De Cruz Date: Fri, 27 Mar 2026 13:46:08 +0100 Subject: [PATCH 5/6] Fix formatting in blending/steps.py --- pysteps/blending/steps.py | 30 +++++++++++++++++------------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index baeff660..2a8aac73 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -1207,9 +1207,10 @@ def transform_to_lagrangian(precip, i): # not strictly necessary but to be consistent with checking only # latest observation rain field to set zero_precip_radar if self.__params.zero_precip_radar: - self.__precip = np.where(np.isfinite(self.__precip), - np.ones(self.__precip.shape)*np.nanmin(self.__precip), - self.__precip + self.__precip = np.where( + np.isfinite(self.__precip), + np.ones(self.__precip.shape) * np.nanmin(self.__precip), + self.__precip, ) # The norain fraction threshold used for nwp is the default value of 0.0, @@ -1490,7 +1491,8 @@ def __estimate_ar_parameters_radar(self): # Finally, transpose GAMMA to ensure that the shape is the same as np.empty((n_cascade_levels, ar_order)) GAMMA = GAMMA.transpose() - if len(GAMMA.shape)==1: GAMMA = GAMMA.reshape((GAMMA.size,1)) + if len(GAMMA.shape) == 1: + GAMMA = GAMMA.reshape((GAMMA.size, 1)) assert GAMMA.shape == ( self.__config.n_cascade_levels, self.__config.ar_order, @@ -3638,7 +3640,7 @@ def forecast( # Check the input precip and ar_order to be consistent # zero-precip in previous time steps has to be removed # (constant field causes autoregression to fail) - precip, ar_order = check_previous_radar_obs(precip,ar_order) + precip, ar_order = check_previous_radar_obs(precip, ar_order) blending_config = StepsBlendingConfig( n_ens_members=n_ens_members, @@ -3702,9 +3704,9 @@ def forecast( # TODO: Where does this piece of code best fit: in utils or inside the class? -def check_previous_radar_obs(precip,ar_order): +def check_previous_radar_obs(precip, ar_order): """Check all radar time steps and remove zero precipitation time steps - + Parameters ---------- precip : array-like @@ -3717,7 +3719,7 @@ def check_previous_radar_obs(precip,ar_order): Returns ------- precip : numpy array - Array of shape (ar_order+1,m,n) containing the modified array with + Array of shape (ar_order+1,m,n) containing the modified array with input precipitation fields ordered by timestamp from oldest to newest. The time steps between the inputs are assumed to be regular. ar_order : int @@ -3728,7 +3730,7 @@ def check_previous_radar_obs(precip,ar_order): assert precip.shape[0] >= 2 # Check all time steps for zero-precip (constant field, minimum==maximum) - zero_precip = [np.nanmin(obs)==np.nanmax(obs) for obs in precip] + zero_precip = [np.nanmin(obs) == np.nanmax(obs) for obs in precip] if zero_precip[-1] or ~np.any(zero_precip): # Unchanged if no rain in latest time step -> will be processed as zero_precip_radar=True # or Unchanged if all time steps contain rain @@ -3739,12 +3741,14 @@ def check_previous_radar_obs(precip,ar_order): # find latest non-zero precip # ATTENTION: This changes the time between precip[-2] and precip[-1] from initial 5min to a longer period prev = np.arange(len(zero_precip[:-2]))[~np.array(zero_precip[:-2])][-1] - return precip[[prev,-1]], 1 - raise ValueError("Precipitation in latest but no previous time step. Not possible to calculate autoregression.") + return precip[[prev, -1]], 1 + raise ValueError( + "Precipitation in latest but no previous time step. Not possible to calculate autoregression." + ) else: # Keep the latest time steps that do all contain precip - precip = precip[np.max(np.arange(len(zero_precip))[zero_precip])+1:].copy() - return precip, min(ar_order,precip.shape[0]-1) + precip = precip[np.max(np.arange(len(zero_precip))[zero_precip]) + 1 :].copy() + return precip, min(ar_order, precip.shape[0] - 1) # TODO: Where does this piece of code best fit: in utils or inside the class? From 4650d8b11610dabd2d511a4af4f0c1df8a218b7b Mon Sep 17 00:00:00 2001 From: FelixE91 Date: Tue, 31 Mar 2026 17:10:00 +0200 Subject: [PATCH 6/6] fix ar_order=1 runs (AR-1 model) --- pysteps/blending/skill_scores.py | 26 +++++++++++++++++++------- pysteps/blending/steps.py | 18 ++++++++++++++---- 2 files changed, 33 insertions(+), 11 deletions(-) diff --git a/pysteps/blending/skill_scores.py b/pysteps/blending/skill_scores.py index c9521de0..c5575f92 100644 --- a/pysteps/blending/skill_scores.py +++ b/pysteps/blending/skill_scores.py @@ -136,7 +136,7 @@ def lt_dependent_cor_nwp(lt, correlations, outdir_path, n_model=0, skill_kwargs= return rho -def lt_dependent_cor_extrapolation(PHI, correlations=None, correlations_prev=None): +def lt_dependent_cor_extrapolation(PHI, correlations=None, correlations_prev=None, ar_order=2): """Determine the correlation of the extrapolation (nowcast) component for lead time lt and cascade k, by assuming that the correlation determined at t=0 regresses towards the climatological values. @@ -153,6 +153,8 @@ def lt_dependent_cor_extrapolation(PHI, correlations=None, correlations_prev=Non be found from the AR-2 model. correlations_prev : array-like, optional Similar to correlations, but from the timestep before that. + ar_order : int, optional + The order of the autoregressive model to use. Must be >= 1. Returns ------- @@ -170,12 +172,22 @@ def lt_dependent_cor_extrapolation(PHI, correlations=None, correlations_prev=Non if correlations_prev is None: correlations_prev = np.repeat(1.0, PHI.shape[0]) # Same for correlations at first time step, we set it to - # phi1 / (1 - phi2), see BPS2004 - if correlations is None: - correlations = PHI[:, 0] / (1.0 - PHI[:, 1]) - - # Calculate the correlation for lead time lt - rho = PHI[:, 0] * correlations + PHI[:, 1] * correlations_prev + if ar_order == 1: + # Yule-Walker equation for AR-1 model -> just PHI + if correlations is None: + correlations = PHI[:, 0] + # Calculate the correlation for lead time lt (AR-1) + rho = PHI[:, 0] * correlations + elif ar_order == 2: + # phi1 / (1 - phi2), see BPS2004 + if correlations is None: + correlations = PHI[:, 0] / (1.0 - PHI[:, 1]) + # Calculate the correlation for lead time lt (AR-2) + rho = PHI[:, 0] * correlations + PHI[:, 1] * correlations_prev + else: + raise ValueError( + "Autoregression of higher order than 2 is not defined. Please set ar_order = 2." + ) # Finally, set the current correlations array as the previous one for the # next time step diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 2a8aac73..89b3bd17 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -1641,13 +1641,23 @@ def __prepare_forecast_loop(self): ) # initizalize the current and previous extrapolation forecast scale for the nowcasting component - # phi1 / (1 - phi2), see BPS2004 + # ar_order=1: rho = phi1 + # ar_order=2: rho = phi1 / (1 - phi2) + # -> see Hamilton1984, Pulkinen2019, BPS2004 (Yule-Walker equations) 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] - ) + if self.__config.ar_order == 1: + self.__state.rho_extrap_cascade = self.__params.PHI[:, 0] + elif self.__config.ar_order == 2: + self.__state.rho_extrap_cascade = self.__params.PHI[:, 0] / ( + 1.0 - self.__params.PHI[:, 1] + ) + else: + raise ValueError( + "Autoregression of higher order than 2 is not defined. Please set ar_order = 2." + ) + def __initialize_noise_cascades(self): """