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'] diff --git a/pysteps/blending/skill_scores.py b/pysteps/blending/skill_scores.py index c9521de0..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): +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 +155,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 +174,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 161c4454..561112a1 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -1196,12 +1196,22 @@ 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 +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)) assert GAMMA.shape == ( self.__config.n_cascade_levels, self.__config.ar_order, @@ -1629,13 +1641,22 @@ 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 Hamilton1994, Pulkkinen2019, 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): """ @@ -2037,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): @@ -2228,9 +2250,12 @@ 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 @@ -3625,6 +3650,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 +3716,65 @@ 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 + 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] + 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. Input delta time changed." + ) + 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() + if precip.shape[0] - 1 < ar_order: + # Give a warning + print( + 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) + + # 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. diff --git a/pysteps/nowcasts/steps.py b/pysteps/nowcasts/steps.py index dc77c7e5..f503aedb 100644 --- a/pysteps/nowcasts/steps.py +++ b/pysteps/nowcasts/steps.py @@ -358,11 +358,19 @@ 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 +1503,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 +1547,61 @@ 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). + """ + 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] + 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. Input delta time changed." + ) + 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() + 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/tests/test_blending_steps.py b/pysteps/tests/test_blending_steps.py index fe9e2863..abfa5522 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,27 +430,127 @@ def test_steps_blending( measure_time=False, ) - assert precip_forecast.ndim == 4, "Wrong amount of dimensions in forecast output" + ### + # The blending + ### + # Test with full radar data + run_and_assert_forecast( + radar_precip, + forecast_kwargs, + expected_n_ens_members, + n_timesteps, + converter, + metadata, + ) - 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" +@pytest.mark.parametrize("ar_order", [1, 2]) +def test_steps_blending_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") - # Transform the data back into mm/h - precip_forecast, _ = converter(precip_forecast, metadata) + 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 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, 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 + 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"] - assert ( - precip_forecast.shape[1] == n_timesteps - ), "Wrong amount of output time steps in converted forecast output" + # Decompose NWP + n_cascade_levels = 6 + 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=1, + 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=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 7e558db4..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 io, motion, nowcasts, verification +from pysteps import io, motion, nowcasts, utils, verification from pysteps.tests.helpers import get_precipitation_fields steps_arg_names = ( @@ -192,3 +192,253 @@ 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( + n_ens_members, + n_cascade_levels, + ar_order, + mask_method, + probmatching_method, + domain, + timesteps, + max_crps, +): + pytest.importorskip("cv2") + + ### + # The input data + ### + # Initialise dummy NWP data + if not isinstance(timesteps, int): + n_timesteps = len(timesteps) + else: + n_timesteps = timesteps + + vel_pert_method = None + + # Define dummy nowcast input data + radar_precip = np.zeros((3, 200, 200)) + + 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" + 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 + + 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" + + ### + # 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, + 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=ar_order, + vel_pert_method=vel_pert_method, + conditional=False, + probmatching_method=probmatching_method, + mask_method=mask_method, + callback=None, + return_output=True, + seed=None, + num_workers=1, + fft_method="numpy", + domain=domain, + 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, + 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, + 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", + callback=None, + return_output=True, + seed=42, + num_workers=1, + fft_method="numpy", + domain="spatial", + 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, + ) 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]