From 84bb003fe9256a9f5893c331d1758c54a83f8451 Mon Sep 17 00:00:00 2001 From: Simon De Kock Date: Wed, 17 Jul 2024 16:48:45 +0200 Subject: [PATCH] Fix of a blending bug due to the removal of rain outside the lagrangian blended mask. Fixed by aplying the smooth mask to the lagrangian blended reference nowcast (R_pm_blended) --- pysteps/blending/steps.py | 42 +++++++++++++++++++++++---------------- pysteps/blending/utils.py | 7 +++++-- 2 files changed, 30 insertions(+), 19 deletions(-) diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index 562298228..020066710 100644 --- a/pysteps/blending/steps.py +++ b/pysteps/blending/steps.py @@ -90,7 +90,7 @@ def forecast( conditional=False, probmatching_method="cdf", mask_method="incremental", - create_smooth_radar_mask=False, + smooth_radar_mask_range=0, callback=None, return_output=True, seed=None, @@ -211,9 +211,10 @@ def forecast( 'obs' = apply precip_thr to the most recently observed precipitation intensity field, 'incremental' = iteratively buffer the mask with a certain rate (currently it is 1 km/min), None=no masking. - create_smooth_radar_mask: bool - If set to True, this generates a smooth mask for the radar that removes - hard edges between NWP and radar in the final blended product. + smooth_radar_mask_range: int, Default is 0. + If 0 this generates a normal forecast. To create a smooth mask, this range + should be a positive value. The smooth radar mask for the radar removes + the hard edges between NWP and radar in the final blended product. probmatching_method: {'cdf','mean',None}, optional Method for matching the statistics of the forecast field with those of the most recently observed one. 'cdf'=map the forecast CDF to the observed @@ -1455,12 +1456,13 @@ def worker(j): # forecast outside the radar domain. Therefore, fill these # areas with the "..._mod_only" blended forecasts, consisting # of the NWP and noise components. - if create_smooth_radar_mask is True: + + nan_indices = np.isnan(R_f_new) + if smooth_radar_mask_range != 0: # Compute the smooth dilated mask - nan_indices = np.isnan(R_f_new) new_mask = blending.utils.compute_smooth_dilated_mask( nan_indices, - 100, + max_padding_size_in_px=smooth_radar_mask_range, gaussian_kernel_size=9, inverted=False, non_linear_growth_kernel_sizes=False, @@ -1484,22 +1486,28 @@ def worker(j): ], axis=0, ) - # Finally, fill the remaining nan values, if present, with - # the minimum value in the forecast - R_f_new[np.isnan(R_f_new)] = np.nanmin(R_f_new) + + nan_indices = np.isnan(R_pm_blended) + R_pm_blended = np.nansum( + [ + R_pm_blended * mask_radar, + R_pm_blended_mod_only * mask_model, + ], + axis=0, + ) else: - nan_indices = np.isnan(R_f_new) R_f_new[nan_indices] = R_f_new_mod_only[nan_indices] nan_indices = np.isnan(R_pm_blended) R_pm_blended[nan_indices] = R_pm_blended_mod_only[ nan_indices ] - # Finally, fill the remaining nan values, if present, with - # the minimum value in the forecast - nan_indices = np.isnan(R_f_new) - R_f_new[nan_indices] = np.nanmin(R_f_new) - nan_indices = np.isnan(R_pm_blended) - R_pm_blended[nan_indices] = np.nanmin(R_pm_blended) + + # Finally, fill the remaining nan values, if present, with + # the minimum value in the forecast + nan_indices = np.isnan(R_f_new) + R_f_new[nan_indices] = np.nanmin(R_f_new) + nan_indices = np.isnan(R_pm_blended) + R_pm_blended[nan_indices] = np.nanmin(R_pm_blended) # 8.7.2. Apply the masking and prob. matching if mask_method is not None: diff --git a/pysteps/blending/utils.py b/pysteps/blending/utils.py index b619e49d0..460e2cdbe 100644 --- a/pysteps/blending/utils.py +++ b/pysteps/blending/utils.py @@ -569,7 +569,7 @@ def check_norain(precip_arr, precip_thr=None, norain_thr=0.0): def compute_smooth_dilated_mask( original_mask, - max_padding_size_in_px, + max_padding_size_in_px=0, gaussian_kernel_size=9, inverted=False, non_linear_growth_kernel_sizes=False, @@ -582,7 +582,7 @@ def compute_smooth_dilated_mask( original_mask : array_like Two-dimensional boolean array containing the input mask. max_padding_size_in_px : int - The maximum size of the padding in pixels. + The maximum size of the padding in pixels. Default is 100. gaussian_kernel_size : int, optional Size of the Gaussian kernel to use for blurring. Default is 9. inverted : bool, optional @@ -603,6 +603,9 @@ def compute_smooth_dilated_mask( " Please install it using `pip install opencv-python`." ) + if max_padding_size_in_px < 0: + raise ValueError("max_padding_size_in_px must be greater than or equal to 0.") + # Convert the original mask to uint8 numpy array and invert if needed array_2d = np.array(original_mask, dtype=np.uint8) if inverted: