From cbd0d06679e41de127ea1100c57050d521e27123 Mon Sep 17 00:00:00 2001 From: FelixE91 Date: Thu, 26 Mar 2026 18:01:25 +0100 Subject: [PATCH 01/20] 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 02/20] 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 03/20] 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 04/20] 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 05/20] 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 06/20] 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): """ From bd7f399d9ee70594e3683c8c11eb30e1996f214d Mon Sep 17 00:00:00 2001 From: FelixE91 Date: Wed, 1 Apr 2026 12:06:48 +0200 Subject: [PATCH 07/20] add warning message for ar_order 1 --- pysteps/blending/steps.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 89b3bd17..754bd831 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -1643,7 +1643,7 @@ def __prepare_forecast_loop(self): # initizalize the current and previous extrapolation forecast scale for the nowcasting component # ar_order=1: rho = phi1 # ar_order=2: rho = phi1 / (1 - phi2) - # -> see Hamilton1984, Pulkinen2019, BPS2004 (Yule-Walker equations) + # -> see Hamilton1994, Pulkkinen2019, BPS2004 (Yule-Walker equations) self.__state.rho_extrap_cascade_prev = np.repeat( 1.0, self.__params.PHI.shape[0] ) @@ -3750,6 +3750,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 + print("[WARNING] Radar input time steps adapted and ar_order set to 1.") prev = np.arange(len(zero_precip[:-2]))[~np.array(zero_precip[:-2])][-1] return precip[[prev, -1]], 1 raise ValueError( From b49902e8dddc2869c6cf2e29ce3ef0befa622f99 Mon Sep 17 00:00:00 2001 From: FelixE91 Date: Wed, 1 Apr 2026 12:27:15 +0200 Subject: [PATCH 08/20] Fix precip-noprecip input array nowcast failure --- pysteps/nowcasts/steps.py | 60 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 59 insertions(+), 1 deletion(-) diff --git a/pysteps/nowcasts/steps.py b/pysteps/nowcasts/steps.py index dc77c7e5..3fb853bd 100644 --- a/pysteps/nowcasts/steps.py +++ b/pysteps/nowcasts/steps.py @@ -358,11 +358,15 @@ def compute_forecast(self): self.__precip = self.__precip[-(self.__config.ar_order + 1) :, :, :].copy() self.__initialize_nowcast_components() if check_norain( - self.__precip, + self.__precip[-1], self.__config.precip_threshold, self.__config.norain_threshold, self.__params.noise_kwargs["win_fun"], ): + # Set all to inputs to 0 (also previous times) as we just check latest input in check_norain + self.__precip = self.__precip.copy() + self.__precip = np.where(np.isfinite(self.__precip),np.ones(self.__precip.shape)*np.nanmin(self.__precip),self.__precip) + return zero_precipitation_forecast( self.__config.n_ens_members, self.__time_steps, @@ -1495,6 +1499,11 @@ def forecast( :cite:`Seed2003`, :cite:`BPS2006`, :cite:`SPN2013`, :cite:`PCH2019b` """ + # 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) + nowcaster_config = StepsNowcasterConfig( n_ens_members=n_ens_members, n_cascade_levels=n_cascade_levels, @@ -1534,3 +1543,52 @@ def forecast( nowcaster.reset_states_and_params() # Call the appropriate methods within the class 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 + print("[WARNING] Radar input time steps adapted and ar_order set to 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: + # 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) From 9135618b22a059f7e49c0c68b65be17a69c04025 Mon Sep 17 00:00:00 2001 From: FelixE91 Date: Wed, 1 Apr 2026 13:06:11 +0200 Subject: [PATCH 09/20] fix mixed precip blending test --- pysteps/tests/test_blending_steps.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/pysteps/tests/test_blending_steps.py b/pysteps/tests/test_blending_steps.py index b2712c68..687284dc 100644 --- a/pysteps/tests/test_blending_steps.py +++ b/pysteps/tests/test_blending_steps.py @@ -446,7 +446,7 @@ def test_steps_blending( @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 + """Test that a forecast succeeds when only the 2 latest radar frames have precipitation (initiating cell corner case that produces NaN autocorrelations for the earlier, all-zero cascade levels).""" pytest.importorskip("cv2") @@ -473,10 +473,14 @@ def test_steps_blending_partial_zero_radar(ar_order): # 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 + radar_precip[-2, 40:125, 30] = 0.5 + radar_precip[-2, 40:125, 31] = 4.5 + radar_precip[-2, 40:125, 32] = 4.0 + radar_precip[-2, 40:125, 33] = 2.0 + radar_precip[-1, 30:155, 32] = 1.0 + radar_precip[-1, 30:155, 33] = 5.0 + radar_precip[-1, 30:155, 34] = 5.0 + radar_precip[-1, 30:155, 35] = 4.5 # Threshold, convert and transform radar_precip[radar_precip < metadata["threshold"]] = 0.0 From 0f45439dd39ae21b8f5ae71c3221f70e36b9926a Mon Sep 17 00:00:00 2001 From: FelixE91 Date: Wed, 1 Apr 2026 13:09:48 +0200 Subject: [PATCH 10/20] Add new nowcast test for mixed precip-no-precip input - similar to the tests in test_blending_steps.py --- pysteps/tests/test_nowcasts_steps.py | 346 ++++++++++++++++++++++++++- 1 file changed, 344 insertions(+), 2 deletions(-) diff --git a/pysteps/tests/test_nowcasts_steps.py b/pysteps/tests/test_nowcasts_steps.py index 7e558db4..42588a33 100644 --- a/pysteps/tests/test_nowcasts_steps.py +++ b/pysteps/tests/test_nowcasts_steps.py @@ -1,10 +1,10 @@ import os -from datetime import timedelta +from datetime import datetime, timedelta import numpy as np import pytest -from pysteps import io, motion, nowcasts, verification +from pysteps import cascade, io, motion, nowcasts, utils, verification from pysteps.tests.helpers import get_precipitation_fields steps_arg_names = ( @@ -192,3 +192,345 @@ def callback(array): timestamps = [startdate + (i + 1) * td for i in range(n_timesteps)] assert (metadata_netcdf["leadtimes"] == leadtimes).all(), "Wrong leadtimes" assert (metadata_netcdf["timestamps"] == timestamps).all(), "Wrong timestamps" + + +def run_and_assert_forecast( + precip, forecast_kwargs, expected_n_ens_members, n_timesteps, converter, metadata +): + """Run a pysteps nowcast and assert the output has the expected shape.""" + precip_forecast = nowcasts.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" + + +@pytest.mark.parametrize(steps_arg_names, steps_arg_values) +def test_steps_nowcast( + timesteps, + n_ens_members, + n_cascade_levels, + nowcasting_method, + mask_method, + probmatching_method, + expected_n_ens_members, + zero_radar, + smooth_radar_mask_range, + resample_distribution, + vel_pert_method, + max_mask_rim, +): + pytest.importorskip("cv2") + + ### + # The input data + ### + # Initialise dummy NWP data + if not isinstance(timesteps, int): + n_timesteps = len(timesteps) + last_timestep = timesteps[-1] + else: + n_timesteps = timesteps + last_timestep = timesteps + + # Define dummy nowcast input data + radar_precip = np.zeros((3, 200, 200)) + + if not zero_radar: + for i in range(2): + radar_precip[i, 5:150, 30 + 1 * i] = 0.1 + radar_precip[i, 5:150, 31 + 1 * i] = 0.5 + radar_precip[i, 5:150, 32 + 1 * i] = 0.5 + radar_precip[i, 5:150, 33 + 1 * i] = 5.0 + radar_precip[i, 5:150, 34 + 1 * i] = 5.0 + radar_precip[i, 5:150, 35 + 1 * i] = 4.5 + radar_precip[i, 5:150, 36 + 1 * i] = 4.5 + radar_precip[i, 5:150, 37 + 1 * i] = 4.0 + radar_precip[i, 5:150, 38 + 1 * i] = 1.0 + radar_precip[i, 5:150, 39 + 1 * i] = 0.5 + radar_precip[i, 5:150, 40 + 1 * i] = 0.5 + radar_precip[i, 5:150, 41 + 1 * i] = 0.1 + radar_precip[2, 30:155, 30 + 1 * 2] = 0.1 + radar_precip[2, 30:155, 31 + 1 * 2] = 0.1 + radar_precip[2, 30:155, 32 + 1 * 2] = 1.0 + radar_precip[2, 30:155, 33 + 1 * 2] = 5.0 + radar_precip[2, 30:155, 34 + 1 * 2] = 5.0 + radar_precip[2, 30:155, 35 + 1 * 2] = 4.5 + radar_precip[2, 30:155, 36 + 1 * 2] = 4.5 + radar_precip[2, 30:155, 37 + 1 * 2] = 4.0 + radar_precip[2, 30:155, 38 + 1 * 2] = 2.0 + radar_precip[2, 30:155, 39 + 1 * 2] = 1.0 + radar_precip[2, 30:155, 40 + 1 * 3] = 0.5 + radar_precip[2, 30:155, 41 + 1 * 3] = 0.1 + + precip_nowcast = np.zeros((n_ens_members, last_timestep + 1, 200, 200)) + + if nowcasting_method == "external_nowcast_ens": + nowcasting_method = "external_nowcast" + for n_ens_member in range(n_ens_members): + for i in range(precip_nowcast.shape[1]): + precip_nowcast[ + n_ens_member, i, 30:165, 30 + 1 * (i + 1) * n_ens_member + ] = 0.1 + precip_nowcast[ + n_ens_member, i, 30:165, 31 + 1 * (i + 1) * n_ens_member + ] = 0.5 + precip_nowcast[ + n_ens_member, i, 30:165, 32 + 1 * (i + 1) * n_ens_member + ] = 0.5 + precip_nowcast[ + n_ens_member, i, 30:165, 33 + 1 * (i + 1) * n_ens_member + ] = 5.0 + precip_nowcast[ + n_ens_member, i, 30:165, 34 + 1 * (i + 1) * n_ens_member + ] = 5.0 + precip_nowcast[ + n_ens_member, i, 30:165, 35 + 1 * (i + 1) * n_ens_member + ] = 4.5 + precip_nowcast[ + n_ens_member, i, 30:165, 36 + 1 * (i + 1) * n_ens_member + ] = 4.5 + precip_nowcast[ + n_ens_member, i, 30:165, 37 + 1 * (i + 1) * n_ens_member + ] = 4.0 + precip_nowcast[ + n_ens_member, i, 30:165, 38 + 1 * (i + 1) * n_ens_member + ] = 1.0 + precip_nowcast[ + n_ens_member, i, 30:165, 39 + 1 * (i + 1) * n_ens_member + ] = 0.5 + precip_nowcast[ + n_ens_member, i, 30:165, 40 + 1 * (i + 1) * n_ens_member + ] = 0.5 + 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" + for i in range(precip_nowcast.shape[1]): + precip_nowcast[0, i, 30:165, 30 + 1 * i] = 0.1 + precip_nowcast[0, i, 30:165, 31 + 1 * i] = 0.5 + precip_nowcast[0, i, 30:165, 32 + 1 * i] = 0.5 + precip_nowcast[0, i, 30:165, 33 + 1 * i] = 5.0 + precip_nowcast[0, i, 30:165, 34 + 1 * i] = 5.0 + precip_nowcast[0, i, 30:165, 35 + 1 * i] = 4.5 + precip_nowcast[0, i, 30:165, 36 + 1 * i] = 4.5 + precip_nowcast[0, i, 30:165, 37 + 1 * i] = 4.0 + precip_nowcast[0, i, 30:165, 38 + 1 * i] = 1.0 + precip_nowcast[0, i, 30:165, 39 + 1 * i] = 0.5 + precip_nowcast[0, i, 30:165, 40 + 1 * i] = 0.5 + precip_nowcast[0, i, 30:165, 41 + 1 * i] = 0.1 + + metadata = dict() + metadata["unit"] = "mm" + metadata["transformation"] = "dB" + metadata["accutime"] = 5.0 + metadata["transform"] = "dB" + metadata["zerovalue"] = 0.0 + metadata["threshold"] = 0.01 + metadata["zr_a"] = 200.0 + metadata["zr_b"] = 1.6 + + # Also set the outdir_path, clim_kwargs and mask_kwargs + outdir_path_skill = "./tmp/" + + if max_mask_rim is not None: + mask_kwargs = dict({"mask_rim": 10, "max_mask_rim": max_mask_rim}) + else: + mask_kwargs = None + + ### + # First threshold the data and convert it to dBR + ### + # threshold the data + radar_precip[radar_precip < metadata["threshold"]] = 0.0 + + # convert the data + converter = utils.get_method("mm/h") + radar_precip, _ = converter(radar_precip, metadata) + + # transform the data + transformer = utils.get_method(metadata["transformation"]) + radar_precip, _ = transformer(radar_precip, metadata) + + # set NaN equal to zero + radar_precip[~np.isfinite(radar_precip)] = metadata["zerovalue"] + + assert ( + np.any(~np.isfinite(radar_precip)) == False + ), "There are still infinite values in the input radar data" + + ### + # Decompose the R_NWP data + ### + + # Initial decomposition settings + decomp_method, _ = cascade.get_method("fft") + bandpass_filter_method = "gaussian" + precip_shape = radar_precip.shape[1:] + filter_method = cascade.get_method(bandpass_filter_method) + bp_filter = filter_method(precip_shape, n_cascade_levels) + + ### + # Determine the velocity fields + ### + oflow_method = motion.get_method("lucaskanade") + radar_velocity = oflow_method(radar_precip) + + ### + # Shared forecast kwargs + ### + forecast_kwargs = dict( + velocity=radar_velocity, + timesteps=timesteps, + timestep=5.0, + issuetime=datetime.strptime("202112012355", "%Y%m%d%H%M"), + n_ens_members=n_ens_members, + n_cascade_levels=n_cascade_levels, + 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, + 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, + mask_kwargs=mask_kwargs, + measure_time=False, + ) + + ### + # The Nowcast + ### + # Test with full radar data + run_and_assert_forecast( + radar_precip, + forecast_kwargs, + expected_n_ens_members, + n_timesteps, + converter, + metadata, + ) + + +@pytest.mark.parametrize("ar_order", [1, 2]) +def test_steps_nowcast_partial_zero_radar(ar_order): + """Test that a forecast succeeds when only the 2 latest radar frames have + precipitation (initiating cell corner case that produces NaN autocorrelations + for the earlier, all-zero cascade levels).""" + pytest.importorskip("cv2") + + 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, + ) + + # Build radar data: only the 2 latest (most recent) frames have precipitation + radar_precip = np.zeros((3, 200, 200)) + radar_precip[-2, 40:125, 30] = 0.5 + radar_precip[-2, 40:125, 31] = 4.5 + radar_precip[-2, 40:125, 32] = 4.0 + radar_precip[-2, 40:125, 33] = 2.0 + radar_precip[-1, 30:155, 32] = 1.0 + radar_precip[-1, 30:155, 33] = 5.0 + radar_precip[-1, 30:155, 34] = 5.0 + radar_precip[-1, 30:155, 35] = 4.5 + + # Threshold, convert and transform + radar_precip[radar_precip < metadata["threshold"]] = 0.0 + converter = utils.get_method("mm/h") + radar_precip, _ = converter(radar_precip, metadata) + transformer = utils.get_method(metadata["transformation"]) + radar_precip, _ = transformer(radar_precip, metadata) + radar_precip[~np.isfinite(radar_precip)] = metadata["zerovalue"] + + + # Velocity fields + oflow_method = motion.get_method("lucaskanade") + radar_velocity = oflow_method(radar_precip) + + run_and_assert_forecast( + radar_precip, + dict( + velocity=radar_velocity, + timesteps=n_timesteps, + timestep=5.0, + issuetime=datetime.strptime("202112012355", "%Y%m%d%H%M"), + n_ens_members=2, + 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, + 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, + mask_kwargs=None, + measure_time=False, + ), + expected_n_ens_members=2, + n_timesteps=n_timesteps, + converter=converter, + metadata=metadata, + ) From 3c723ce27f3fc66941e4f2250a22b7d263e05144 Mon Sep 17 00:00:00 2001 From: FelixE91 Date: Wed, 1 Apr 2026 14:35:29 +0200 Subject: [PATCH 11/20] improve warning messages in check_previous_radar_obs --- pysteps/blending/steps.py | 5 ++++- pysteps/nowcasts/steps.py | 5 ++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 754bd831..06024d69 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -3750,7 +3750,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 - print("[WARNING] Radar input time steps adapted and ar_order set to 1.") + print("[WARNING] Radar input time steps adapted and ar_order set to 1. Input delta time changed.") prev = np.arange(len(zero_precip[:-2]))[~np.array(zero_precip[:-2])][-1] return precip[[prev, -1]], 1 raise ValueError( @@ -3759,6 +3759,9 @@ def check_previous_radar_obs(precip, ar_order): else: # Keep the latest time steps that do all contain precip precip = precip[np.max(np.arange(len(zero_precip))[zero_precip]) + 1 :].copy() + if precip.shape[0]-1 < ar_order: + # Give a warning + print(f"[WARNING] Radar input only with {precip.shape[0]} non-zero time steps and ar_order set to {precip.shape[0]-1}.") return precip, min(ar_order, precip.shape[0] - 1) diff --git a/pysteps/nowcasts/steps.py b/pysteps/nowcasts/steps.py index 3fb853bd..20cc2c54 100644 --- a/pysteps/nowcasts/steps.py +++ b/pysteps/nowcasts/steps.py @@ -1582,7 +1582,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 - print("[WARNING] Radar input time steps adapted and ar_order set to 1.") + print("[WARNING] Radar input time steps adapted and ar_order set to 1. Input delta time changed.") prev = np.arange(len(zero_precip[:-2]))[~np.array(zero_precip[:-2])][-1] return precip[[prev, -1]], 1 raise ValueError( @@ -1591,4 +1591,7 @@ def check_previous_radar_obs(precip, ar_order): else: # Keep the latest time steps that do all contain precip precip = precip[np.max(np.arange(len(zero_precip))[zero_precip]) + 1 :].copy() + if precip.shape[0]-1 < ar_order: + # Give a warning + print(f"[WARNING] Radar input only with {precip.shape[0]} non-zero time steps and ar_order set to {precip.shape[0]-1}.") return precip, min(ar_order, precip.shape[0] - 1) From 7d3a2df0ce0516ecd1e65e12e53a85d080214e1b Mon Sep 17 00:00:00 2001 From: FelixE91 Date: Wed, 1 Apr 2026 16:29:07 +0200 Subject: [PATCH 12/20] fix renormalize cascade --- pysteps/blending/steps.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 06024d69..e90f90a4 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -2250,9 +2250,8 @@ def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): ) ) # Renormalize the cascade - worker_state.precip_cascades[j][i][1] /= np.std( - worker_state.precip_cascades[j][i][1] - ) + for nl, obs_cascade_level in enumerate(worker_state.precip_cascades[j][i]): + worker_state.precip_cascades[j][i][nl] /= np.std(obs_cascade_level) else: # use the deterministic AR(p) model computed above if # perturbations are disabled From 2cb883f41afc3841d8e8ef421bc19f767d5c5b9e Mon Sep 17 00:00:00 2001 From: Lesley De Cruz Date: Wed, 1 Apr 2026 16:28:18 +0200 Subject: [PATCH 13/20] Fix n_ens_members to 1 in blending tests and update nowcast test parameters inherited from blending test --- pysteps/tests/test_blending_steps.py | 4 ++-- pysteps/tests/test_nowcasts_steps.py | 31 ++++++++++++++-------------- 2 files changed, 17 insertions(+), 18 deletions(-) diff --git a/pysteps/tests/test_blending_steps.py b/pysteps/tests/test_blending_steps.py index 687284dc..c3daf7c5 100644 --- a/pysteps/tests/test_blending_steps.py +++ b/pysteps/tests/test_blending_steps.py @@ -518,7 +518,7 @@ def test_steps_blending_partial_zero_radar(ar_order): timesteps=n_timesteps, timestep=5.0, issuetime=datetime.datetime.strptime("202112012355", "%Y%m%d%H%M"), - n_ens_members=2, + n_ens_members=1, n_cascade_levels=n_cascade_levels, blend_nwp_members=False, precip_thr=metadata["threshold"], @@ -551,7 +551,7 @@ def test_steps_blending_partial_zero_radar(ar_order): mask_kwargs=None, measure_time=False, ), - expected_n_ens_members=2, + expected_n_ens_members=1, n_timesteps=n_timesteps, converter=converter, metadata=metadata, diff --git a/pysteps/tests/test_nowcasts_steps.py b/pysteps/tests/test_nowcasts_steps.py index 42588a33..c7276ed0 100644 --- a/pysteps/tests/test_nowcasts_steps.py +++ b/pysteps/tests/test_nowcasts_steps.py @@ -1,5 +1,5 @@ import os -from datetime import datetime, timedelta +from datetime import timedelta import numpy as np import pytest @@ -224,18 +224,14 @@ def run_and_assert_forecast( @pytest.mark.parametrize(steps_arg_names, steps_arg_values) def test_steps_nowcast( - timesteps, n_ens_members, n_cascade_levels, - nowcasting_method, + ar_order, mask_method, probmatching_method, - expected_n_ens_members, - zero_radar, - smooth_radar_mask_range, - resample_distribution, - vel_pert_method, - max_mask_rim, + domain, + timesteps, + max_crps, ): pytest.importorskip("cv2") @@ -250,6 +246,14 @@ def test_steps_nowcast( n_timesteps = timesteps last_timestep = timesteps + nowcasting_method = "steps" + expected_n_ens_members = n_ens_members + zero_radar = False + smooth_radar_mask_range = 0 + resample_distribution = False + vel_pert_method = None + max_mask_rim = None + # Define dummy nowcast input data radar_precip = np.zeros((3, 200, 200)) @@ -404,7 +408,6 @@ def test_steps_nowcast( velocity=radar_velocity, timesteps=timesteps, timestep=5.0, - issuetime=datetime.strptime("202112012355", "%Y%m%d%H%M"), n_ens_members=n_ens_members, n_cascade_levels=n_cascade_levels, precip_thr=metadata["threshold"], @@ -414,7 +417,7 @@ def test_steps_nowcast( bandpass_filter_method="gaussian", noise_method="nonparametric", noise_stddev_adj="auto", - ar_order=2, + ar_order=ar_order, vel_pert_method=vel_pert_method, conditional=False, probmatching_method=probmatching_method, @@ -426,8 +429,7 @@ def test_steps_nowcast( seed=None, num_workers=1, fft_method="numpy", - domain="spatial", - outdir_path_skill=outdir_path_skill, + domain=domain, extrap_kwargs=None, filter_kwargs=None, noise_kwargs=None, @@ -488,7 +490,6 @@ def test_steps_nowcast_partial_zero_radar(ar_order): radar_precip, _ = transformer(radar_precip, metadata) radar_precip[~np.isfinite(radar_precip)] = metadata["zerovalue"] - # Velocity fields oflow_method = motion.get_method("lucaskanade") radar_velocity = oflow_method(radar_precip) @@ -499,7 +500,6 @@ def test_steps_nowcast_partial_zero_radar(ar_order): velocity=radar_velocity, timesteps=n_timesteps, timestep=5.0, - issuetime=datetime.strptime("202112012355", "%Y%m%d%H%M"), n_ens_members=2, precip_thr=metadata["threshold"], kmperpixel=1.0, @@ -521,7 +521,6 @@ def test_steps_nowcast_partial_zero_radar(ar_order): num_workers=1, fft_method="numpy", domain="spatial", - outdir_path_skill="./tmp/", extrap_kwargs=None, filter_kwargs=None, noise_kwargs=None, From 26fa133d89bfa2e335c3e2028046bd21b504ce7a Mon Sep 17 00:00:00 2001 From: Lesley De Cruz Date: Wed, 1 Apr 2026 16:45:06 +0200 Subject: [PATCH 14/20] Remove more leftovers from blending from the test_nowcasts_steps.py script. --- pysteps/tests/test_nowcasts_steps.py | 147 +++++---------------------- 1 file changed, 28 insertions(+), 119 deletions(-) diff --git a/pysteps/tests/test_nowcasts_steps.py b/pysteps/tests/test_nowcasts_steps.py index c7276ed0..c1ef73dd 100644 --- a/pysteps/tests/test_nowcasts_steps.py +++ b/pysteps/tests/test_nowcasts_steps.py @@ -4,7 +4,7 @@ import numpy as np import pytest -from pysteps import cascade, io, motion, nowcasts, utils, verification +from pysteps import io, motion, nowcasts, utils, verification from pysteps.tests.helpers import get_precipitation_fields steps_arg_names = ( @@ -241,109 +241,39 @@ def test_steps_nowcast( # Initialise dummy NWP data if not isinstance(timesteps, int): n_timesteps = len(timesteps) - last_timestep = timesteps[-1] else: n_timesteps = timesteps - last_timestep = timesteps - nowcasting_method = "steps" - expected_n_ens_members = n_ens_members - zero_radar = False - smooth_radar_mask_range = 0 - resample_distribution = False vel_pert_method = None - max_mask_rim = None # Define dummy nowcast input data radar_precip = np.zeros((3, 200, 200)) - if not zero_radar: - for i in range(2): - radar_precip[i, 5:150, 30 + 1 * i] = 0.1 - radar_precip[i, 5:150, 31 + 1 * i] = 0.5 - radar_precip[i, 5:150, 32 + 1 * i] = 0.5 - radar_precip[i, 5:150, 33 + 1 * i] = 5.0 - radar_precip[i, 5:150, 34 + 1 * i] = 5.0 - radar_precip[i, 5:150, 35 + 1 * i] = 4.5 - radar_precip[i, 5:150, 36 + 1 * i] = 4.5 - radar_precip[i, 5:150, 37 + 1 * i] = 4.0 - radar_precip[i, 5:150, 38 + 1 * i] = 1.0 - radar_precip[i, 5:150, 39 + 1 * i] = 0.5 - radar_precip[i, 5:150, 40 + 1 * i] = 0.5 - radar_precip[i, 5:150, 41 + 1 * i] = 0.1 - radar_precip[2, 30:155, 30 + 1 * 2] = 0.1 - radar_precip[2, 30:155, 31 + 1 * 2] = 0.1 - radar_precip[2, 30:155, 32 + 1 * 2] = 1.0 - radar_precip[2, 30:155, 33 + 1 * 2] = 5.0 - radar_precip[2, 30:155, 34 + 1 * 2] = 5.0 - radar_precip[2, 30:155, 35 + 1 * 2] = 4.5 - radar_precip[2, 30:155, 36 + 1 * 2] = 4.5 - radar_precip[2, 30:155, 37 + 1 * 2] = 4.0 - radar_precip[2, 30:155, 38 + 1 * 2] = 2.0 - radar_precip[2, 30:155, 39 + 1 * 2] = 1.0 - radar_precip[2, 30:155, 40 + 1 * 3] = 0.5 - radar_precip[2, 30:155, 41 + 1 * 3] = 0.1 - - precip_nowcast = np.zeros((n_ens_members, last_timestep + 1, 200, 200)) - - if nowcasting_method == "external_nowcast_ens": - nowcasting_method = "external_nowcast" - for n_ens_member in range(n_ens_members): - for i in range(precip_nowcast.shape[1]): - precip_nowcast[ - n_ens_member, i, 30:165, 30 + 1 * (i + 1) * n_ens_member - ] = 0.1 - precip_nowcast[ - n_ens_member, i, 30:165, 31 + 1 * (i + 1) * n_ens_member - ] = 0.5 - precip_nowcast[ - n_ens_member, i, 30:165, 32 + 1 * (i + 1) * n_ens_member - ] = 0.5 - precip_nowcast[ - n_ens_member, i, 30:165, 33 + 1 * (i + 1) * n_ens_member - ] = 5.0 - precip_nowcast[ - n_ens_member, i, 30:165, 34 + 1 * (i + 1) * n_ens_member - ] = 5.0 - precip_nowcast[ - n_ens_member, i, 30:165, 35 + 1 * (i + 1) * n_ens_member - ] = 4.5 - precip_nowcast[ - n_ens_member, i, 30:165, 36 + 1 * (i + 1) * n_ens_member - ] = 4.5 - precip_nowcast[ - n_ens_member, i, 30:165, 37 + 1 * (i + 1) * n_ens_member - ] = 4.0 - precip_nowcast[ - n_ens_member, i, 30:165, 38 + 1 * (i + 1) * n_ens_member - ] = 1.0 - precip_nowcast[ - n_ens_member, i, 30:165, 39 + 1 * (i + 1) * n_ens_member - ] = 0.5 - precip_nowcast[ - n_ens_member, i, 30:165, 40 + 1 * (i + 1) * n_ens_member - ] = 0.5 - 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" - for i in range(precip_nowcast.shape[1]): - precip_nowcast[0, i, 30:165, 30 + 1 * i] = 0.1 - precip_nowcast[0, i, 30:165, 31 + 1 * i] = 0.5 - precip_nowcast[0, i, 30:165, 32 + 1 * i] = 0.5 - precip_nowcast[0, i, 30:165, 33 + 1 * i] = 5.0 - precip_nowcast[0, i, 30:165, 34 + 1 * i] = 5.0 - precip_nowcast[0, i, 30:165, 35 + 1 * i] = 4.5 - precip_nowcast[0, i, 30:165, 36 + 1 * i] = 4.5 - precip_nowcast[0, i, 30:165, 37 + 1 * i] = 4.0 - precip_nowcast[0, i, 30:165, 38 + 1 * i] = 1.0 - precip_nowcast[0, i, 30:165, 39 + 1 * i] = 0.5 - precip_nowcast[0, i, 30:165, 40 + 1 * i] = 0.5 - precip_nowcast[0, i, 30:165, 41 + 1 * i] = 0.1 + for i in range(2): + radar_precip[i, 5:150, 30 + 1 * i] = 0.1 + radar_precip[i, 5:150, 31 + 1 * i] = 0.5 + radar_precip[i, 5:150, 32 + 1 * i] = 0.5 + radar_precip[i, 5:150, 33 + 1 * i] = 5.0 + radar_precip[i, 5:150, 34 + 1 * i] = 5.0 + radar_precip[i, 5:150, 35 + 1 * i] = 4.5 + radar_precip[i, 5:150, 36 + 1 * i] = 4.5 + radar_precip[i, 5:150, 37 + 1 * i] = 4.0 + radar_precip[i, 5:150, 38 + 1 * i] = 1.0 + radar_precip[i, 5:150, 39 + 1 * i] = 0.5 + radar_precip[i, 5:150, 40 + 1 * i] = 0.5 + radar_precip[i, 5:150, 41 + 1 * i] = 0.1 + radar_precip[2, 30:155, 30 + 1 * 2] = 0.1 + radar_precip[2, 30:155, 31 + 1 * 2] = 0.1 + radar_precip[2, 30:155, 32 + 1 * 2] = 1.0 + radar_precip[2, 30:155, 33 + 1 * 2] = 5.0 + radar_precip[2, 30:155, 34 + 1 * 2] = 5.0 + radar_precip[2, 30:155, 35 + 1 * 2] = 4.5 + radar_precip[2, 30:155, 36 + 1 * 2] = 4.5 + radar_precip[2, 30:155, 37 + 1 * 2] = 4.0 + radar_precip[2, 30:155, 38 + 1 * 2] = 2.0 + radar_precip[2, 30:155, 39 + 1 * 2] = 1.0 + radar_precip[2, 30:155, 40 + 1 * 3] = 0.5 + radar_precip[2, 30:155, 41 + 1 * 3] = 0.1 metadata = dict() metadata["unit"] = "mm" @@ -355,13 +285,7 @@ def test_steps_nowcast( metadata["zr_a"] = 200.0 metadata["zr_b"] = 1.6 - # Also set the outdir_path, clim_kwargs and mask_kwargs - outdir_path_skill = "./tmp/" - - if max_mask_rim is not None: - mask_kwargs = dict({"mask_rim": 10, "max_mask_rim": max_mask_rim}) - else: - mask_kwargs = None + mask_kwargs = None ### # First threshold the data and convert it to dBR @@ -384,17 +308,6 @@ def test_steps_nowcast( np.any(~np.isfinite(radar_precip)) == False ), "There are still infinite values in the input radar data" - ### - # Decompose the R_NWP data - ### - - # Initial decomposition settings - decomp_method, _ = cascade.get_method("fft") - bandpass_filter_method = "gaussian" - precip_shape = radar_precip.shape[1:] - filter_method = cascade.get_method(bandpass_filter_method) - bp_filter = filter_method(precip_shape, n_cascade_levels) - ### # Determine the velocity fields ### @@ -422,8 +335,6 @@ def test_steps_nowcast( 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, @@ -445,7 +356,7 @@ def test_steps_nowcast( run_and_assert_forecast( radar_precip, forecast_kwargs, - expected_n_ens_members, + n_ens_members, n_timesteps, converter, metadata, @@ -513,8 +424,6 @@ def test_steps_nowcast_partial_zero_radar(ar_order): conditional=False, probmatching_method=None, mask_method="incremental", - resample_distribution=False, - smooth_radar_mask_range=0, callback=None, return_output=True, seed=42, From e605150a9bdf8cd769e5c0df6c9fdb2487529c17 Mon Sep 17 00:00:00 2001 From: FelixE91 Date: Thu, 2 Apr 2026 18:11:56 +0200 Subject: [PATCH 15/20] reformatted to match default --- pysteps/blending/steps.py | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index e90f90a4..6d776fbc 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -1656,8 +1656,7 @@ def __prepare_forecast_loop(self): else: raise ValueError( "Autoregression of higher order than 2 is not defined. Please set ar_order = 2." - ) - + ) def __initialize_noise_cascades(self): """ @@ -2059,6 +2058,7 @@ def __determine_skill_for_current_timestep(self, t): PHI=self.__params.PHI, correlations=self.__state.rho_extrap_cascade, correlations_prev=self.__state.rho_extrap_cascade_prev, + ar_order=self.__config.ar_order, ) def __determine_NWP_skill_for_next_timestep(self, t, j, worker_state): @@ -2250,8 +2250,12 @@ def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t): ) ) # Renormalize the cascade - for nl, obs_cascade_level in enumerate(worker_state.precip_cascades[j][i]): - worker_state.precip_cascades[j][i][nl] /= np.std(obs_cascade_level) + for nl, obs_cascade_level in enumerate( + worker_state.precip_cascades[j][i] + ): + worker_state.precip_cascades[j][i][nl] /= np.std( + obs_cascade_level + ) else: # use the deterministic AR(p) model computed above if # perturbations are disabled @@ -3749,7 +3753,9 @@ 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 - print("[WARNING] Radar input time steps adapted and ar_order set to 1. Input delta time changed.") + print( + "[WARNING] Radar input time steps adapted and ar_order set to 1. Input delta time changed." + ) prev = np.arange(len(zero_precip[:-2]))[~np.array(zero_precip[:-2])][-1] return precip[[prev, -1]], 1 raise ValueError( @@ -3758,9 +3764,11 @@ def check_previous_radar_obs(precip, ar_order): else: # Keep the latest time steps that do all contain precip precip = precip[np.max(np.arange(len(zero_precip))[zero_precip]) + 1 :].copy() - if precip.shape[0]-1 < ar_order: + if precip.shape[0] - 1 < ar_order: # Give a warning - print(f"[WARNING] Radar input only with {precip.shape[0]} non-zero time steps and ar_order set to {precip.shape[0]-1}.") + print( + f"[WARNING] Radar input only with {precip.shape[0]} non-zero time steps and ar_order set to {precip.shape[0]-1}." + ) return precip, min(ar_order, precip.shape[0] - 1) From c47d49d420b7c4a0c0d80429b0e0feca24a3a0ce Mon Sep 17 00:00:00 2001 From: FelixE91 Date: Fri, 3 Apr 2026 09:18:15 +0200 Subject: [PATCH 16/20] black reformatting --- pysteps/blending/skill_scores.py | 4 +++- pysteps/nowcasts/steps.py | 16 ++++++++++++---- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/pysteps/blending/skill_scores.py b/pysteps/blending/skill_scores.py index c5575f92..8bdccaae 100644 --- a/pysteps/blending/skill_scores.py +++ b/pysteps/blending/skill_scores.py @@ -136,7 +136,9 @@ 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, ar_order=2): +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. diff --git a/pysteps/nowcasts/steps.py b/pysteps/nowcasts/steps.py index 20cc2c54..f0160c88 100644 --- a/pysteps/nowcasts/steps.py +++ b/pysteps/nowcasts/steps.py @@ -365,7 +365,11 @@ def compute_forecast(self): ): # Set all to inputs to 0 (also previous times) as we just check latest input in check_norain self.__precip = self.__precip.copy() - 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, + ) return zero_precipitation_forecast( self.__config.n_ens_members, @@ -1582,7 +1586,9 @@ 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 - print("[WARNING] Radar input time steps adapted and ar_order set to 1. Input delta time changed.") + print( + "[WARNING] Radar input time steps adapted and ar_order set to 1. Input delta time changed." + ) prev = np.arange(len(zero_precip[:-2]))[~np.array(zero_precip[:-2])][-1] return precip[[prev, -1]], 1 raise ValueError( @@ -1591,7 +1597,9 @@ def check_previous_radar_obs(precip, ar_order): else: # Keep the latest time steps that do all contain precip precip = precip[np.max(np.arange(len(zero_precip))[zero_precip]) + 1 :].copy() - if precip.shape[0]-1 < ar_order: + if precip.shape[0] - 1 < ar_order: # Give a warning - print(f"[WARNING] Radar input only with {precip.shape[0]} non-zero time steps and ar_order set to {precip.shape[0]-1}.") + print( + f"[WARNING] Radar input only with {precip.shape[0]} non-zero time steps and ar_order set to {precip.shape[0]-1}." + ) return precip, min(ar_order, precip.shape[0] - 1) From 3188d38876a4371a3108c7158dbc30f0df4d0a4c Mon Sep 17 00:00:00 2001 From: FelixE91 Date: Fri, 3 Apr 2026 12:55:46 +0200 Subject: [PATCH 17/20] black reformatted --- pysteps/blending/steps.py | 7 +++++-- pysteps/tests/test_blending_steps.py | 2 -- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 6d776fbc..561112a1 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -3740,7 +3740,10 @@ def check_previous_radar_obs(precip, ar_order): 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 + if not precip.shape[0] >= 2: + raise ValueError( + "Wrong precip shape. The radar input must have at least 2 time steps." + ) # Check all time steps for zero-precip (constant field, minimum==maximum) zero_precip = [np.nanmin(obs) == np.nanmax(obs) for obs in precip] @@ -3767,7 +3770,7 @@ def check_previous_radar_obs(precip, ar_order): if precip.shape[0] - 1 < ar_order: # Give a warning print( - f"[WARNING] Radar input only with {precip.shape[0]} non-zero time steps and ar_order set to {precip.shape[0]-1}." + f"[WARNING] Remove zero-precip time steps from radar input.\nRemaining non-zero time steps: {precip.shape[0]}, ar_order set to {precip.shape[0]-1}." ) return precip, min(ar_order, precip.shape[0] - 1) diff --git a/pysteps/tests/test_blending_steps.py b/pysteps/tests/test_blending_steps.py index c3daf7c5..abfa5522 100644 --- a/pysteps/tests/test_blending_steps.py +++ b/pysteps/tests/test_blending_steps.py @@ -496,8 +496,6 @@ def test_steps_blending_partial_zero_radar(ar_order): # 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 From 1e05808d5ba3bffdd32a9676234a1e3dbd1e99c3 Mon Sep 17 00:00:00 2001 From: FelixE91 Date: Fri, 3 Apr 2026 13:07:59 +0200 Subject: [PATCH 18/20] fix codacy issue --- pysteps/nowcasts/steps.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pysteps/nowcasts/steps.py b/pysteps/nowcasts/steps.py index f0160c88..f503aedb 100644 --- a/pysteps/nowcasts/steps.py +++ b/pysteps/nowcasts/steps.py @@ -1572,8 +1572,10 @@ def check_previous_radar_obs(precip, ar_order): 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 + if not precip.shape[0] >= 2: + raise ValueError( + "Wrong precip shape. The radar input must have at least 2 time steps." + ) # Check all time steps for zero-precip (constant field, minimum==maximum) zero_precip = [np.nanmin(obs) == np.nanmax(obs) for obs in precip] From 2cc64a4d9c2492d6ce43dc87ced40d9bbb5fd6b8 Mon Sep 17 00:00:00 2001 From: FelixE91 Date: Fri, 3 Apr 2026 14:37:01 +0200 Subject: [PATCH 19/20] add bandit config to exclude assert check in tests --- .pre-commit-config.yaml | 7 +++++++ bandit.yaml | 3 +++ 2 files changed, 10 insertions(+) create mode 100644 bandit.yaml diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 10cd6dbc..7050fa3c 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -4,3 +4,10 @@ repos: hooks: - id: black language_version: python3 + +- repo: https://github.com/PyCQA/bandit + rev: 1.9.4 + hooks: + - id: bandit + args: ["-c", "bandit.yaml"] + additional_dependencies: ["bandit[yaml]"] diff --git a/bandit.yaml b/bandit.yaml new file mode 100644 index 00000000..cc15b0f9 --- /dev/null +++ b/bandit.yaml @@ -0,0 +1,3 @@ +# skip check for assert usage in the tests +assert_used: + skips: ['*test_*.py'] From 12d10a5f8c711efe984c4539bedacff8b4a8ac0d Mon Sep 17 00:00:00 2001 From: FelixE91 Date: Fri, 3 Apr 2026 17:45:05 +0200 Subject: [PATCH 20/20] Improve code readability --- pysteps/timeseries/autoregression.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/pysteps/timeseries/autoregression.py b/pysteps/timeseries/autoregression.py index 0617fd3a..7ce75e81 100644 --- a/pysteps/timeseries/autoregression.py +++ b/pysteps/timeseries/autoregression.py @@ -425,7 +425,7 @@ def estimate_ar_params_yw(gamma, d=0, check_stationarity=True): Returns ------- out: ndarray - Array of length p+1 containing the AR(p) parameters for for the + Array of length p+1 containing the AR(p) parameters for the lag-p terms and the innovation term. Notes @@ -455,15 +455,16 @@ def estimate_ar_params_yw(gamma, d=0, check_stationarity=True): "Error in estimate_ar_params_yw: " "nonstationary AR(p) process" ) + # Estimate PHI_{j,0} as phi_pert. See Pulkkinen et al. (2019), eq. 6. c = 1.0 for j in range(p): c -= gamma[j] * phi[j] - phi_pert = np.sqrt(c) - # If the expression inside the square root is negative, phi_pert cannot # be computed and it is set to zero instead. - if not np.isfinite(phi_pert): + if c < 0: phi_pert = 0.0 + else: + phi_pert = np.sqrt(c) if d == 1: phi = _compute_differenced_model_params(phi, p, 1, 1) @@ -499,7 +500,7 @@ def estimate_ar_params_yw_localized(gamma, d=0): Returns ------- out: list - List of length p+1 containing the AR(p) parameter fields for for the + List of length p+1 containing the AR(p) parameter fields for the lag-p terms and the innovation term. The parameter fields have the same shape as the elements of gamma. @@ -538,6 +539,7 @@ def estimate_ar_params_yw_localized(gamma, d=0): phi[:, i] = phi_ + # Estimate PHI_{j,0} as phi_pert. See Pulkkinen et al. (2019), eq. 6. c = 1.0 for i in range(p): c -= gamma_1d[i] * phi[i]