Skip to content

Use max rain pixels for timestep and member to decompose for precip cascades in case of zero radar #461

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Mar 31, 2025
Merged

Use max rain pixels for timestep and member to decompose for precip cascades in case of zero radar #461

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
41b621a
Use max rain pixels for timestep and member to decompose for precip c…
mats-knmi Feb 25, 2025
2a0d956
store precip threshold in state
mats-knmi Mar 3, 2025
3cbd296
Merge branch 'master' into zero-radar-precip-cascades-max-pixels
dnerini Mar 28, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
92 changes: 41 additions & 51 deletions pysteps/blending/steps.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -325,6 +325,7 @@
precip_models_provided_is_cascade: bool = False
xy_coordinates: np.ndarray | None = None
precip_zerovalue: float | None = None
precip_threshold: float | None = None
mask_threshold: np.ndarray | None = None
zero_precip_radar: bool = False
zero_precip_model_fields: bool = False
Expand Down Expand Up @@ -806,6 +807,8 @@
else:
self.__params.mask_kwargs = deepcopy(self.__config.mask_kwargs)

self.__params.precip_threshold = self.__config.precip_threshold

if np.any(~np.isfinite(self.__velocity)):
raise ValueError("velocity contains non-finite values")

Expand All @@ -815,12 +818,12 @@
% self.__config.mask_method
)

if self.__config.conditional and self.__config.precip_threshold is None:
if self.__config.conditional and self.__params.precip_threshold is None:
raise ValueError("conditional=True but precip_thr is not set")

if (
self.__config.mask_method is not None
and self.__config.precip_threshold is None
and self.__params.precip_threshold is None
):
raise ValueError("mask_method!=None but precip_thr=None")

Expand Down Expand Up @@ -931,7 +934,7 @@
) = (None, None)

if self.__config.conditional or self.__config.mask_method is not None:
print(f"precip. intensity threshold: {self.__config.precip_threshold}")
print(f"precip. intensity threshold: {self.__params.precip_threshold}")
print(f"no-rain fraction threshold for radar: {self.__config.norain_threshold}")
print("")

Expand Down Expand Up @@ -994,7 +997,7 @@
# TODO: is this logical_and correct here? Now only those places where precip is in all images is saved?
self.__params.mask_threshold = np.logical_and.reduce(
[
self.__precip[i, :, :] >= self.__config.precip_threshold
self.__precip[i, :, :] >= self.__params.precip_threshold
for i in range(self.__precip.shape[0])
]
)
Expand Down Expand Up @@ -1098,15 +1101,15 @@
# Check for zero input fields in the radar and NWP data.
self.__params.zero_precip_radar = check_norain(
self.__precip,
self.__config.precip_threshold,
self.__params.precip_threshold,
self.__config.norain_threshold,
self.__params.noise_kwargs["win_fun"],
)
# The norain fraction threshold used for nwp is the default value of 0.0,
# since nwp does not suffer from clutter.
self.__params.zero_precip_model_fields = check_norain(
self.__precip_models,
self.__config.precip_threshold,
self.__params.precip_threshold,
self.__config.norain_threshold,
self.__params.noise_kwargs["win_fun"],
)
Expand Down Expand Up @@ -1173,43 +1176,6 @@
updates the cascade with NWP data, uses the NWP velocity field, and
initializes the noise model based on the time step with the most rain.
"""
# If zero_precip_radar, make sure that precip_cascade does not contain
# only nans or infs. If so, fill it with the zero value.

# Look for a timestep and member with rain so that we have a sensible decomposition
done = False
for t in self.__timesteps:
if done:
break
for j in range(self.__precip_models.shape[0]):
if not check_norain(
self.__precip_models[j, t],
self.__config.precip_threshold,
self.__config.norain_threshold,
self.__params.noise_kwargs["win_fun"],
):
if self.__state.precip_models_cascades is not None:
self.__state.precip_cascades[
~np.isfinite(self.__state.precip_cascades)
] = np.nanmin(
self.__state.precip_models_cascades[j, t]["cascade_levels"]
)
continue
precip_models_cascade_timestep = self.__params.decomposition_method(
self.__precip_models[j, t, :, :],
bp_filter=self.__params.bandpass_filter,
fft_method=self.__params.fft,
output_domain=self.__config.domain,
normalize=True,
compute_stats=True,
compact_output=True,
)["cascade_levels"]
self.__state.precip_cascades[
~np.isfinite(self.__state.precip_cascades)
] = np.nanmin(precip_models_cascade_timestep)
done = True
break

# If zero_precip_radar is True, only use the velocity field of the NWP
# forecast. I.e., velocity (radar) equals velocity_model at the first time
# step.
Expand All @@ -1227,16 +1193,16 @@
# might be zero as well). Else, initialize the noise with the radar
# rainfall data
# Initialize noise based on the NWP field time step where the fraction of rainy cells is highest
if self.__config.precip_threshold is None:
self.__config.precip_threshold = np.nanmin(self.__precip_models)
if self.__params.precip_threshold is None:
self.__params.precip_threshold = np.nanmin(self.__precip_models)

Check warning on line 1197 in pysteps/blending/steps.py

View check run for this annotation

Codecov / codecov/patch

pysteps/blending/steps.py#L1197 

Added line #L1197 was not covered by tests

max_rain_pixels = -1
max_rain_pixels_j = -1
max_rain_pixels_t = -1
for j in range(self.__precip_models.shape[0]):
for t in self.__timesteps:
rain_pixels = self.__precip_models[j][t][
self.__precip_models[j][t] > self.__config.precip_threshold
self.__precip_models[j][t] > self.__params.precip_threshold
].size
if rain_pixels > max_rain_pixels:
max_rain_pixels = rain_pixels
Expand All @@ -1249,6 +1215,30 @@
np.float64, copy=False
)

# If zero_precip_radar, make sure that precip_cascade does not contain
# only nans or infs. If so, fill it with the zero value.
if self.__state.precip_models_cascades is not None:
self.__state.precip_cascades[~np.isfinite(self.__state.precip_cascades)] = (
np.nanmin(
self.__state.precip_models_cascades[
max_rain_pixels_j, max_rain_pixels_t
]["cascade_levels"]
)
)
else:
precip_models_cascade_timestep = self.__params.decomposition_method(
self.__precip_models[max_rain_pixels_j, max_rain_pixels_t, :, :],
bp_filter=self.__params.bandpass_filter,
fft_method=self.__params.fft,
output_domain=self.__config.domain,
normalize=True,
compute_stats=True,
compact_output=True,
)["cascade_levels"]
self.__state.precip_cascades[~np.isfinite(self.__state.precip_cascades)] = (
np.nanmin(precip_models_cascade_timestep)
)

# Make sure precip_noise_input is three-dimensional
if len(self.__state.precip_noise_input.shape) != 3:
self.__state.precip_noise_input = self.__state.precip_noise_input[
Expand Down Expand Up @@ -1281,7 +1271,7 @@
precip_forecast_min = np.min(self.__state.precip_noise_input)
self.__params.noise_std_coeffs = noise.utils.compute_noise_stddev_adjs(
self.__state.precip_noise_input[-1, :, :],
self.__config.precip_threshold,
self.__params.precip_threshold,
precip_forecast_min,
self.__params.bandpass_filter,
self.__params.decomposition_method,
Expand Down Expand Up @@ -2687,7 +2677,7 @@
# steps, the buffer mask keeps increasing.
precip_field_mask = (
precip_forecast_probability_matching_blended
>= self.__config.precip_threshold
>= self.__params.precip_threshold
)

# Buffer the mask
Expand Down Expand Up @@ -2731,7 +2721,7 @@
# rainfall field
precip_field_mask_temp = (
precip_forecast_probability_matching_blended
>= self.__config.precip_threshold
>= self.__params.precip_threshold
)

# Set to min value outside of mask
Expand Down Expand Up @@ -2789,12 +2779,12 @@
mean_probabiltity_matching_forecast = np.mean(
precip_forecast_probability_matching_resampled[
precip_forecast_probability_matching_resampled
>= self.__config.precip_threshold
>= self.__params.precip_threshold
]
)
no_rain_mask = (
worker_state.final_blended_forecast_recomposed
>= self.__config.precip_threshold
>= self.__params.precip_threshold
)
mean_precip_forecast = np.mean(
worker_state.final_blended_forecast_recomposed[no_rain_mask]
Expand Down