Allow for semi dynamic incremental mask in blending setup

pysteps/blending/steps.py

@@ -147,10 +147,10 @@ class StepsBlendingConfig:
       field, 'incremental' = iteratively buffer the mask with a certain rate
       (currently it is 1 km/min), None=no masking.
     resample_distribution: bool, optional
-        Method to resample the distribution from the extrapolation and NWP cascade as input
-        for the probability matching. Not resampling these distributions may lead to losing
-        some extremes when the weight of both the extrapolation and NWP cascade is similar.
-        Defaults to True.
+      Method to resample the distribution from the extrapolation and NWP cascade as input
+      for the probability matching. Not resampling these distributions may lead to losing
+      some extremes when the weight of both the extrapolation and NWP cascade is similar.
+      Defaults to True.
     smooth_radar_mask_range: int, Default is 0.
       Method to smooth the transition between the radar-NWP-noise blend and the NWP-noise
       blend near the edge of the radar domain (radar mask), where the radar data is either
@@ -249,9 +249,9 @@ class StepsBlendingConfig:
       (the number of NWP models) and 'window_length' (the minimum number of
       days the clim file should have, otherwise the default is used).
     mask_kwargs: dict
-      Optional dictionary containing mask keyword arguments 'mask_f' and
-      'mask_rim', the factor defining the the mask increment and the rim size,
-      respectively.
+      Optional dictionary containing mask keyword arguments 'mask_f',
+      'mask_rim' and 'max_mask_rim', the factor defining the the mask
+      increment and the (maximum) rim size, respectively.
       The mask increment is defined as mask_f*timestep/kmperpixel.
     measure_time: bool
       If set to True, measure, print and return the computation time.
@@ -639,7 +639,10 @@ def worker(j):
                         self.__recompose_cascade_to_rainfall_field(j, worker_state)
                         final_blended_forecast_single_member = (
                             self.__post_process_output(
-                                j, final_blended_forecast_single_member, worker_state
+                                j,
+                                t_sub,
+                                final_blended_forecast_single_member,
+                                worker_state,
                             )
                         )
                     final_blended_forecast_all_members_one_timestep[j] = (
@@ -1480,6 +1483,9 @@ def __prepare_forecast_loop(self):
         if self.__config.mask_method == "incremental":
             # get mask parameters
             self.__params.mask_rim = self.__params.mask_kwargs.get("mask_rim", 10)
+            self.__params.max_mask_rim = self.__params.mask_kwargs.get(
+                "max_mask_rim", 10
+            )
             mask_f = self.__params.mask_kwargs.get("mask_f", 1.0)
             # initialize the structuring element
             struct = generate_binary_structure(2, 1)
@@ -2511,7 +2517,7 @@ def __recompose_cascade_to_rainfall_field(self, j, worker_state):
             )
 
     def __post_process_output(
-        self, j, final_blended_forecast_single_member, worker_state
+        self, j, t_sub, final_blended_forecast_single_member, worker_state
     ):
         """
         Apply post-processing steps to refine the final blended forecast. This
@@ -2674,16 +2680,16 @@ def __post_process_output(
                 worker_state.final_blended_forecast_recomposed.min()
             )
             if self.__config.mask_method == "incremental":
-                # The incremental mask is slightly different from
-                # the implementation in the non-blended steps.py, as
-                # it is not based on the last forecast, but instead
-                # on R_pm_blended. Therefore, the buffer does not
-                # increase over time.
-                # Get the mask for this forecast
+                # The incremental mask is slightly different from the implementation in
+                # nowcasts.steps.py, as it is not computed in the Lagrangian space. Instead,
+                # we use precip_forecast_probability_matched and let the mask_rim increase with
+                # the time step until mask_rim_max. This ensures that for the first t time
+                # steps, the buffer mask keeps increasing.
                 precip_field_mask = (
                     precip_forecast_probability_matching_blended
                     >= self.__config.precip_threshold
                 )
+
                 # Buffer the mask
                 # Convert the precipitation field mask into an 8-bit unsigned integer mask
                 obs_mask_uint8 = precip_field_mask.astype("uint8")
@@ -2698,7 +2704,10 @@ def __post_process_output(
                 accumulated_mask = dilated_mask.astype(float)
 
                 # Iteratively dilate the mask and accumulate the results to create a grayscale rim
-                for _ in range(self.__params.mask_rim):
+                mask_rim_temp = min(
+                    self.__params.mask_rim + t_sub - 1, self.__params.max_mask_rim
+                )
+                for _ in range(mask_rim_temp):
                     dilated_mask = binary_dilation(dilated_mask, struct_element)
                     accumulated_mask += dilated_mask