fix noise tapering and error handling (#460)

* fix noise tapering and error handling

* slight reduction in skill

* add generic check_norain method

* Add test for all-zero steps nowcast (this will fail)

* extend error message

* Add test for no rain to all nowcast methods

* Fix tests

* More fixes for tests

* fix  most nowcasts for norain

* fix linda nowcast

* update error message and docstring

* fix black

* set default check norain win fun to none

* re insert accidentally removed notes from sseps

* Update pysteps/blending/utils.py

Co-authored-by: Daniele Nerini <[email protected]>

---------

Co-authored-by: Daniele Nerini <[email protected]>

Author: mats-knmi

pysteps/blending/steps.py

@@ -56,6 +56,7 @@
 from pysteps.nowcasts import utils as nowcast_utils
 from pysteps.postprocessing import probmatching
 from pysteps.timeseries import autoregression, correlation
+from pysteps.utils.check_norain import check_norain
 
 try:
     import dask
@@ -77,7 +78,7 @@ class StepsBlendingConfig:
     precip_threshold: float, optional
       Specifies the threshold value for minimum observable precipitation
       intensity. Required if mask_method is not None or conditional is True.
-    norain_threshold: float, optional
+    norain_threshold: float
       Specifies the threshold value for the fraction of rainy (see above) pixels
       in the radar rainfall field below which we consider there to be no rain.
       Depends on the amount of clutter typically present.
@@ -435,7 +436,6 @@ def __init__(
 
         # Additional variables for time measurement
         self.__start_time_init = None
-        self.__zero_precip_time = None
         self.__init_time = None
         self.__mainloop_time = None
 
@@ -777,7 +777,7 @@ def __check_inputs(self):
             self.__params.filter_kwargs = deepcopy(self.__config.filter_kwargs)
 
         if self.__config.noise_kwargs is None:
-            self.__params.noise_kwargs = dict()
+            self.__params.noise_kwargs = {"win_fun": "tukey"}
         else:
             self.__params.noise_kwargs = deepcopy(self.__config.noise_kwargs)
 
@@ -1093,17 +1093,19 @@ def transform_to_lagrangian(precip, i):
             self.__precip_models = np.stack(temp_precip_models)
 
         # Check for zero input fields in the radar and NWP data.
-        self.__params.zero_precip_radar = blending.utils.check_norain(
+        self.__params.zero_precip_radar = check_norain(
             self.__precip,
             self.__config.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 = blending.utils.check_norain(
+        self.__params.zero_precip_model_fields = check_norain(
             self.__precip_models,
             self.__config.precip_threshold,
             self.__config.norain_threshold,
+            self.__params.noise_kwargs["win_fun"],
         )
 
     def __zero_precipitation_forecast(self):
@@ -1141,7 +1143,7 @@ def __zero_precipitation_forecast(self):
                 precip_forecast_workers = None
 
         if self.__config.measure_time:
-            self.__zero_precip_time = time.time() - self.__start_time_init
+            zero_precip_time = time.time() - self.__start_time_init
 
         if self.__config.return_output:
             precip_forecast_all_members_all_times = np.stack(
@@ -1154,8 +1156,8 @@ def __zero_precipitation_forecast(self):
             if self.__config.measure_time:
                 return (
                     precip_forecast_all_members_all_times,
-                    self.__zero_precip_time,
-                    self.__zero_precip_time,
+                    zero_precip_time,
+                    zero_precip_time,
                 )
             else:
                 return precip_forecast_all_members_all_times
@@ -1177,10 +1179,11 @@ def __prepare_nowcast_for_zero_radar(self):
             if done:
                 break
             for j in range(self.__precip_models.shape[0]):
-                if not blending.utils.check_norain(
+                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[
@@ -2925,7 +2928,7 @@ def forecast(
     precip_thr: float, optional
       Specifies the threshold value for minimum observable precipitation
       intensity. Required if mask_method is not None or conditional is True.
-    norain_thr: float, optional
+    norain_thr: float
       Specifies the threshold value for the fraction of rainy (see above) pixels
       in the radar rainfall field below which we consider there to be no rain.
       Depends on the amount of clutter typically present.

pysteps/blending/utils.py

@@ -15,11 +15,11 @@
     decompose_NWP
     compute_store_nwp_motion
     load_NWP
-    check_norain
     compute_smooth_dilated_mask
 """
 
 import datetime
+import warnings
 from pathlib import Path
 
 import numpy as np
@@ -28,6 +28,7 @@
 from pysteps.cascade.bandpass_filters import filter_gaussian
 from pysteps.exceptions import MissingOptionalDependency
 from pysteps.utils import get_method as utils_get_method
+from pysteps.utils.check_norain import check_norain as new_check_norain
 
 try:
     import netCDF4
@@ -534,7 +535,7 @@ def load_NWP(input_nc_path_decomp, input_path_velocities, start_time, n_timestep
 
 def check_norain(precip_arr, precip_thr=None, norain_thr=0.0):
     """
-
+    DEPRECATED use :py:mod:`pysteps.utils.check_norain.check_norain` in stead
     Parameters
     ----------
     precip_arr:  array-like
@@ -551,15 +552,10 @@ def check_norain(precip_arr, precip_thr=None, norain_thr=0.0):
       Returns whether the fraction of rainy pixels is below the norain_thr threshold.
 
     """
-
-    if precip_thr is None:
-        precip_thr = np.nanmin(precip_arr)
-    rain_pixels = precip_arr[precip_arr > precip_thr]
-    norain = rain_pixels.size / precip_arr.size <= norain_thr
-    print(
-        f"Rain fraction is: {str(rain_pixels.size / precip_arr.size)}, while minimum fraction is {str(norain_thr)}"
+    warnings.warn(
+        "pysteps.blending.utils.check_norain has been deprecated, use pysteps.utils.check_norain.check_norain instead"
     )
-    return norain
+    return new_check_norain(precip_arr, precip_thr, norain_thr, None)
 
 
 def compute_smooth_dilated_mask(

pysteps/noise/fftgenerators.py

@@ -46,6 +46,7 @@
 
 import numpy as np
 from scipy import optimize
+
 from .. import utils
 
 
@@ -99,7 +100,13 @@ def initialize_param_2d_fft_filter(field, **kwargs):
     if len(field.shape) < 2 or len(field.shape) > 3:
         raise ValueError("the input is not two- or three-dimensional array")
     if np.any(~np.isfinite(field)):
-        raise ValueError("field contains non-finite values")
+        raise ValueError(
+            "field contains non-finite values, this typically happens when the input\n"
+            + "precipitation field provided to pysteps contains (mostly)zero values.\n"
+            + "To prevent this error please call pysteps.utils.check_norain first,\n"
+            + "using the same win_fun as used in this method (tukey by default)\n"
+            + "and then only call this method if that check fails."
+        )
 
     # defaults
     win_fun = kwargs.get("win_fun", None)
@@ -254,7 +261,13 @@ def initialize_nonparam_2d_fft_filter(field, **kwargs):
     if len(field.shape) < 2 or len(field.shape) > 3:
         raise ValueError("the input is not two- or three-dimensional array")
     if np.any(~np.isfinite(field)):
-        raise ValueError("field contains non-finite values")
+        raise ValueError(
+            "field contains non-finite values, this typically happens when the input\n"
+            + "precipitation field provided to pysteps contains (mostly)zero values.\n"
+            + "To prevent this error please call pysteps.utils.check_norain first,\n"
+            + "using the same win_fun as used in this method (tukey by default)\n"
+            + "and then only call this method if that check fails."
+        )
 
     # defaults
     win_fun = kwargs.get("win_fun", "tukey")
@@ -361,7 +374,13 @@ def generate_noise_2d_fft_filter(
     if len(F.shape) != 2:
         raise ValueError("field is not two-dimensional array")
     if np.any(~np.isfinite(F)):
-        raise ValueError("field contains non-finite values")
+        raise ValueError(
+            "field contains non-finite values, this typically happens when the input\n"
+            + "precipitation field provided to pysteps contains (mostly)zero values.\n"
+            + "To prevent this error please call pysteps.utils.check_norain first,\n"
+            + "using the same win_fun as used in this method (tukey by default)\n"
+            + "and then only call this method if that check fails."
+        )
 
     if randstate is None:
         randstate = np.random
@@ -755,7 +774,13 @@ def generate_noise_2d_ssft_filter(F, randstate=None, seed=None, **kwargs):
     if len(F.shape) != 4:
         raise ValueError("the input is not four-dimensional array")
     if np.any(~np.isfinite(F)):
-        raise ValueError("field contains non-finite values")
+        raise ValueError(
+            "field contains non-finite values, this typically happens when the input\n"
+            + "precipitation field provided to pysteps contains (mostly) zero value.s\n"
+            + "To prevent this error please call pysteps.utils.check_norain first,\n"
+            + "using the same win_fun as used in this method (tukey by default)\n"
+            + "and then only call this method if that check fails."
+        )
 
     if "domain" in kwargs.keys() and kwargs["domain"] == "spectral":
         raise NotImplementedError(

pysteps/nowcasts/linda.py

@@ -40,21 +40,23 @@
 import time
 import warnings
 
+from pysteps.utils.check_norain import check_norain
+
 try:
     import dask
 
     DASK_IMPORTED = True
 except ImportError:
     DASK_IMPORTED = False
 import numpy as np
+from scipy import optimize as opt
+from scipy import stats
 from scipy.integrate import nquad
 from scipy.interpolate import interp1d
-from scipy import optimize as opt
 from scipy.signal import convolve
-from scipy import stats
 
 from pysteps import extrapolation, feature, noise
-from pysteps.nowcasts.utils import nowcast_main_loop
+from pysteps.nowcasts.utils import nowcast_main_loop, zero_precipitation_forecast
 
 
 def forecast(
@@ -292,6 +294,19 @@ def forecast(
         True if np.any(~np.isfinite(precip)) else False
     )
 
+    starttime_init = time.time()
+
+    if check_norain(precip, 0.0, 0.0, None):
+        return zero_precipitation_forecast(
+            n_ens_members if nowcast_type == "ensemble" else None,
+            timesteps,
+            precip,
+            callback,
+            return_output,
+            measure_time,
+            starttime_init,
+        )
+
     forecast_gen = _linda_deterministic_init(
         precip,
         velocity,

pysteps/nowcasts/sprog.py

@@ -10,16 +10,16 @@
     forecast
 """
 
-import numpy as np
 import time
 
-from pysteps import cascade
-from pysteps import extrapolation
-from pysteps import utils
+import numpy as np
+
+from pysteps import cascade, extrapolation, utils
 from pysteps.nowcasts import utils as nowcast_utils
+from pysteps.nowcasts.utils import compute_percentile_mask, nowcast_main_loop
 from pysteps.postprocessing import probmatching
 from pysteps.timeseries import autoregression, correlation
-from pysteps.nowcasts.utils import compute_percentile_mask, nowcast_main_loop
+from pysteps.utils.check_norain import check_norain
 
 try:
     import dask
@@ -34,6 +34,7 @@ def forecast(
     velocity,
     timesteps,
     precip_thr=None,
+    norain_thr=0.0,
     n_cascade_levels=6,
     extrap_method="semilagrangian",
     decomp_method="fft",
@@ -68,6 +69,11 @@ def forecast(
         of the list are required to be in ascending order.
     precip_thr: float, required
         The threshold value for minimum observable precipitation intensity.
+    norain_thr: float
+      Specifies the threshold value for the fraction of rainy (see above) pixels
+      in the radar rainfall field below which we consider there to be no rain.
+      Depends on the amount of clutter typically present.
+      Standard set to 0.0
     n_cascade_levels: int, optional
         The number of cascade levels to use. Defaults to 6, see issue #385
         on GitHub.
@@ -182,6 +188,8 @@ def forecast(
 
     if measure_time:
         starttime_init = time.time()
+    else:
+        starttime_init = None
 
     fft = utils.get_method(fft_method, shape=precip.shape[1:], n_threads=num_workers)
 
@@ -203,6 +211,11 @@ def forecast(
         [~np.isfinite(precip[i, :]) for i in range(precip.shape[0])]
     )
 
+    if check_norain(precip, precip_thr, norain_thr, None):
+        return nowcast_utils.zero_precipitation_forecast(
+            None, timesteps, precip, None, True, measure_time, starttime_init
+        )
+
     # determine the precipitation threshold mask
     if conditional:
         mask_thr = np.logical_and.reduce(

pysteps/nowcasts/sseps.py

@@ -18,17 +18,16 @@
     forecast
 """
 
-import numpy as np
 import time
-from scipy.ndimage import generate_binary_structure, iterate_structure
 
+import numpy as np
+from scipy.ndimage import generate_binary_structure, iterate_structure
 
-from pysteps import cascade
-from pysteps import extrapolation
-from pysteps import noise
+from pysteps import cascade, extrapolation, noise
 from pysteps.nowcasts import utils as nowcast_utils
 from pysteps.postprocessing import probmatching
 from pysteps.timeseries import autoregression, correlation
+from pysteps.utils.check_norain import check_norain
 
 try:
     import dask
@@ -211,7 +210,7 @@ def forecast(
         filter_kwargs = dict()
 
     if noise_kwargs is None:
-        noise_kwargs = dict()
+        noise_kwargs = {"win_fun": "tukey"}
 
     if vel_pert_kwargs is None:
         vel_pert_kwargs = dict()
@@ -297,6 +296,8 @@ def forecast(
 
     if measure_time:
         starttime_init = time.time()
+    else:
+        starttime_init = None
 
     # get methods
     extrapolator_method = extrapolation.get_method(extrap_method)
@@ -312,6 +313,22 @@ def forecast(
     if noise_method is not None:
         init_noise, generate_noise = noise.get_method(noise_method)
 
+    if check_norain(
+        precip,
+        precip_thr,
+        war_thr,
+        noise_kwargs["win_fun"],
+    ):
+        return nowcast_utils.zero_precipitation_forecast(
+            n_ens_members,
+            timesteps,
+            precip,
+            callback,
+            return_output,
+            measure_time,
+            starttime_init,
+        )
+
     # advect the previous precipitation fields to the same position with the
     # most recent one (i.e. transform them into the Lagrangian coordinates)
     precip = precip[-(ar_order + 1) :, :, :].copy()