From 4c43aa40a9e9d5765557d62a9b1dead1496378d1 Mon Sep 17 00:00:00 2001
From: mats-knmi <145579783+mats-knmi@users.noreply.github.com>
Date: Mon, 3 Mar 2025 15:43:23 +0100
Subject: [PATCH] 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 <daniele.nerini@gmail.com>

---------

Co-authored-by: Daniele Nerini <daniele.nerini@gmail.com>
---
 pysteps/blending/steps.py            |  23 +++--
 pysteps/blending/utils.py            |  16 ++--
 pysteps/noise/fftgenerators.py       |  33 ++++++-
 pysteps/nowcasts/linda.py            |  21 +++-
 pysteps/nowcasts/sprog.py            |  23 ++++-
 pysteps/nowcasts/sseps.py            |  29 ++++--
 pysteps/nowcasts/steps.py            |  52 +++++++---
 pysteps/nowcasts/utils.py            | 137 ++++++++++++++++++++++++---
 pysteps/tests/test_blending_utils.py |  34 +++----
 pysteps/tests/test_nowcasts_anvil.py |  23 +++++
 pysteps/tests/test_nowcasts_linda.py |  26 +++++
 pysteps/tests/test_nowcasts_sprog.py |  24 +++++
 pysteps/tests/test_nowcasts_sseps.py |  32 +++++++
 pysteps/tests/test_nowcasts_steps.py |  30 +++++-
 pysteps/utils/check_norain.py        |  51 ++++++++++
 pysteps/utils/tapering.py            |  16 ++--
 16 files changed, 483 insertions(+), 87 deletions(-)
 create mode 100644 pysteps/utils/check_norain.py

diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py
index 740fce809..c722e436d 100644
--- a/pysteps/blending/steps.py
+++ b/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.
diff --git a/pysteps/blending/utils.py b/pysteps/blending/utils.py
index 486d7af88..aaed2cfa2 100644
--- a/pysteps/blending/utils.py
+++ b/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(
diff --git a/pysteps/noise/fftgenerators.py b/pysteps/noise/fftgenerators.py
index 7ffabdc26..d3414d06d 100644
--- a/pysteps/noise/fftgenerators.py
+++ b/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(
diff --git a/pysteps/nowcasts/linda.py b/pysteps/nowcasts/linda.py
index 5bd35e1f9..cee990221 100644
--- a/pysteps/nowcasts/linda.py
+++ b/pysteps/nowcasts/linda.py
@@ -40,6 +40,8 @@
 import time
 import warnings
 
+from pysteps.utils.check_norain import check_norain
+
 try:
     import dask
 
@@ -47,14 +49,14 @@
 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,
diff --git a/pysteps/nowcasts/sprog.py b/pysteps/nowcasts/sprog.py
index 3d9e78fa0..3742556e2 100644
--- a/pysteps/nowcasts/sprog.py
+++ b/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(
diff --git a/pysteps/nowcasts/sseps.py b/pysteps/nowcasts/sseps.py
index 94cd34570..f2891b522 100644
--- a/pysteps/nowcasts/sseps.py
+++ b/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()
diff --git a/pysteps/nowcasts/steps.py b/pysteps/nowcasts/steps.py
index 818123da4..dc77c7e59 100644
--- a/pysteps/nowcasts/steps.py
+++ b/pysteps/nowcasts/steps.py
@@ -11,22 +11,24 @@
     forecast
 """
 
-import numpy as np
-from scipy.ndimage import generate_binary_structure, iterate_structure
 import time
 from copy import deepcopy
+from dataclasses import dataclass, field
+from typing import Any, Callable
+
+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 utils
+from pysteps import cascade, extrapolation, noise, utils
 from pysteps.nowcasts import utils as nowcast_utils
+from pysteps.nowcasts.utils import (
+    compute_percentile_mask,
+    nowcast_main_loop,
+    zero_precipitation_forecast,
+)
 from pysteps.postprocessing import probmatching
 from pysteps.timeseries import autoregression, correlation
-from pysteps.nowcasts.utils import compute_percentile_mask, nowcast_main_loop
-
-from dataclasses import dataclass, field
-from typing import Any, Callable
+from pysteps.utils.check_norain import check_norain
 
 try:
     import dask
@@ -50,6 +52,11 @@ class StepsNowcasterConfig:
     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
+      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
     kmperpixel: float, optional
         Spatial resolution of the input data (kilometers/pixel). Required if
         vel_pert_method is not None or mask_method is 'incremental'.
@@ -201,6 +208,7 @@ class StepsNowcasterConfig:
     n_ens_members: int = 24
     n_cascade_levels: int = 6
     precip_threshold: float | None = None
+    norain_threshold: float = 0.0
     kmperpixel: float | None = None
     timestep: float | None = None
     extrapolation_method: str = "semilagrangian"
@@ -349,6 +357,21 @@ def compute_forecast(self):
         # Slice the precipitation field to only use the last ar_order + 1 fields
         self.__precip = self.__precip[-(self.__config.ar_order + 1) :, :, :].copy()
         self.__initialize_nowcast_components()
+        if check_norain(
+            self.__precip,
+            self.__config.precip_threshold,
+            self.__config.norain_threshold,
+            self.__params.noise_kwargs["win_fun"],
+        ):
+            return zero_precipitation_forecast(
+                self.__config.n_ens_members,
+                self.__time_steps,
+                self.__precip,
+                self.__config.callback,
+                self.__config.return_output,
+                self.__config.measure_time,
+                self.__start_time_init,
+            )
 
         self.__perform_extrapolation()
         self.__apply_noise_and_ar_model()
@@ -501,7 +524,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)
 
@@ -1246,6 +1269,7 @@ def forecast(
     n_ens_members=24,
     n_cascade_levels=6,
     precip_thr=None,
+    norain_thr=0.0,
     kmperpixel=None,
     timestep=None,
     extrap_method="semilagrangian",
@@ -1297,6 +1321,11 @@ 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
+      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
     kmperpixel: float, optional
         Spatial resolution of the input data (kilometers/pixel). Required if
         vel_pert_method is not None or mask_method is 'incremental'.
@@ -1470,6 +1499,7 @@ def forecast(
         n_ens_members=n_ens_members,
         n_cascade_levels=n_cascade_levels,
         precip_threshold=precip_thr,
+        norain_threshold=norain_thr,
         kmperpixel=kmperpixel,
         timestep=timestep,
         extrapolation_method=extrap_method,
diff --git a/pysteps/nowcasts/utils.py b/pysteps/nowcasts/utils.py
index fd111e28d..8ddd3da0f 100644
--- a/pysteps/nowcasts/utils.py
+++ b/pysteps/nowcasts/utils.py
@@ -17,6 +17,7 @@
 """
 
 import time
+
 import numpy as np
 from scipy.ndimage import binary_dilation, generate_binary_structure
 
@@ -137,6 +138,130 @@ def compute_percentile_mask(precip, pct):
     return precip >= precip_pct_thr
 
 
+def zero_precipitation_forecast(
+    n_ens_members,
+    timesteps,
+    precip,
+    callback,
+    return_output,
+    measure_time,
+    start_time_init,
+):
+    """
+    Generate a zero-precipitation forecast (filled with the minimum precip value)
+    when no precipitation above the threshold is detected. The forecast is
+    optionally returned or passed to a callback.
+
+    Parameters
+    ----------
+    n_ens_members: int, optional
+        The number of ensemble members to generate.
+    timesteps: int or list of floats
+        Number of time steps to forecast or a list of time steps for which the
+        forecasts are computed (relative to the input time step). The elements
+        of the list are required to be in ascending order.
+    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.
+    callback: function, optional
+        Optional function that is called after computation of each time step of
+        the nowcast. The function takes one argument: a three-dimensional array
+        of shape (n_ens_members,h,w), where h and w are the height and width
+        of the input precipitation fields, respectively. This can be used, for
+        instance, writing the outputs into files.
+    return_output: bool, optional
+        Set to False to disable returning the outputs as numpy arrays. This can
+        save memory if the intermediate results are written to output files using
+        the callback function.
+    measure_time: bool
+        If set to True, measure, print and return the computation time.
+    start_time_init: float
+        The value of the start time counter used to compute total run time.
+
+    Returns
+    -------
+    out: ndarray
+        If return_output is True, a four-dimensional array of shape
+        (n_ens_members,num_timesteps,m,n) containing a time series of forecast
+        precipitation fields for each ensemble member. Otherwise, a None value
+        is returned. The time series starts from t0+timestep, where timestep is
+        taken from the input precipitation fields. If measure_time is True, the
+        return value is a three-element tuple containing the nowcast array, the
+        initialization time of the nowcast generator and the time used in the
+        main loop (seconds).
+    """
+    print("No precipitation above the threshold found in the radar field")
+    print("The resulting forecast will contain only zeros")
+    return_single_member = False
+    if n_ens_members is None:
+        n_ens_members = 1
+        return_single_member = True
+    # Create the output list
+    precip_forecast = [[] for j in range(n_ens_members)]
+
+    # Save per time step to ensure the array does not become too large if
+    # no return_output is requested and callback is not None.
+    timesteps, _, __ = create_timestep_range(timesteps)
+    for t, subtimestep_idx in enumerate(timesteps):
+        # If the timestep is not the first one, we need to provide the zero forecast
+        if t > 0:
+            # Create an empty np array with shape [n_ens_members, rows, cols]
+            # and fill it with the minimum value from precip (corresponding to
+            # zero precipitation)
+            N, M = precip.shape[1:]
+            precip_forecast_workers = np.full((n_ens_members, N, M), np.nanmin(precip))
+            if subtimestep_idx:
+                if callback is not None:
+                    if precip_forecast_workers.shape[1] > 0:
+                        callback(precip_forecast_workers.squeeze())
+                if return_output:
+                    for j in range(n_ens_members):
+                        precip_forecast[j].append(precip_forecast_workers[j])
+            precip_forecast_workers = None
+
+    if measure_time:
+        zero_precip_time = time.time() - start_time_init
+
+    if return_output:
+        precip_forecast_all_members_all_times = np.stack(
+            [np.stack(precip_forecast[j]) for j in range(n_ens_members)]
+        )
+        if return_single_member:
+            precip_forecast_all_members_all_times = (
+                precip_forecast_all_members_all_times[0]
+            )
+
+        if measure_time:
+            return (
+                precip_forecast_all_members_all_times,
+                zero_precip_time,
+                zero_precip_time,
+            )
+        else:
+            return precip_forecast_all_members_all_times
+    else:
+        return None
+
+
+def create_timestep_range(timesteps):
+    """
+    create a range of time steps
+    if an integer time step is given, create a simple range iterator
+    otherwise, assing the time steps to integer bins so that each bin
+    contains a list of time steps belonging to that bin
+    """
+    if isinstance(timesteps, int):
+        timesteps = range(timesteps + 1)
+        timestep_type = "int"
+        original_timesteps = None
+    else:
+        original_timesteps = [0] + list(timesteps)
+        timesteps = binned_timesteps(original_timesteps)
+        timestep_type = "list"
+    return timesteps, original_timesteps, timestep_type
+
+
 def nowcast_main_loop(
     precip,
     velocity,
@@ -219,17 +344,7 @@ def nowcast_main_loop(
     """
     precip_forecast_out = None
 
-    # create a range of time steps
-    # if an integer time step is given, create a simple range iterator
-    # otherwise, assing the time steps to integer bins so that each bin
-    # contains a list of time steps belonging to that bin
-    if isinstance(timesteps, int):
-        timesteps = range(timesteps + 1)
-        timestep_type = "int"
-    else:
-        original_timesteps = [0] + list(timesteps)
-        timesteps = binned_timesteps(original_timesteps)
-        timestep_type = "list"
+    timesteps, original_timesteps, timestep_type = create_timestep_range(timesteps)
 
     state_cur = state
     if not ensemble:
diff --git a/pysteps/tests/test_blending_utils.py b/pysteps/tests/test_blending_utils.py
index 22d0a0045..401b6f1ce 100644
--- a/pysteps/tests/test_blending_utils.py
+++ b/pysteps/tests/test_blending_utils.py
@@ -3,21 +3,21 @@
 import os
 
 import numpy as np
-from numpy.testing import assert_array_almost_equal
 import pytest
+from numpy.testing import assert_array_almost_equal
 
 import pysteps
 from pysteps.blending.utils import (
-    stack_cascades,
     blend_cascades,
-    recompose_cascade,
     blend_optical_flows,
-    decompose_NWP,
+    compute_smooth_dilated_mask,
     compute_store_nwp_motion,
+    decompose_NWP,
     load_NWP,
-    check_norain,
-    compute_smooth_dilated_mask,
+    recompose_cascade,
+    stack_cascades,
 )
+from pysteps.utils.check_norain import check_norain
 
 pytest.importorskip("netCDF4")
 
@@ -398,28 +398,30 @@ def test_blending_utils(
 
     precip_arr = precip_nwp
     # rainy fraction is 0.005847
-    assert not check_norain(precip_arr)
-    assert not check_norain(precip_arr, precip_thr=nwp_metadata["threshold"])
+    assert not check_norain(precip_arr, win_fun=None)
+    assert not check_norain(
+        precip_arr, precip_thr=nwp_metadata["threshold"], win_fun=None
+    )
     assert not check_norain(
-        precip_arr, precip_thr=nwp_metadata["threshold"], norain_thr=0.005
+        precip_arr, precip_thr=nwp_metadata["threshold"], norain_thr=0.005, win_fun=None
     )
-    assert not check_norain(precip_arr, norain_thr=0.005)
+    assert not check_norain(precip_arr, norain_thr=0.005, win_fun=None)
     # so with norain_thr beyond this number it should report that there's no rain
-    assert check_norain(precip_arr, norain_thr=0.006)
+    assert check_norain(precip_arr, norain_thr=0.006, win_fun=None)
     assert check_norain(
-        precip_arr, precip_thr=nwp_metadata["threshold"], norain_thr=0.006
+        precip_arr, precip_thr=nwp_metadata["threshold"], norain_thr=0.006, win_fun=None
     )
 
     # also if we set the precipitation threshold sufficiently high, it should report there's no rain
     # rainy fraction > 4mm/h is 0.004385
-    assert not check_norain(precip_arr, precip_thr=4.0, norain_thr=0.004)
-    assert check_norain(precip_arr, precip_thr=4.0, norain_thr=0.005)
+    assert not check_norain(precip_arr, precip_thr=4.0, norain_thr=0.004, win_fun=None)
+    assert check_norain(precip_arr, precip_thr=4.0, norain_thr=0.005, win_fun=None)
 
     # no rain above 100mm/h so it should give norain
-    assert check_norain(precip_arr, precip_thr=100)
+    assert check_norain(precip_arr, precip_thr=100, win_fun=None)
 
     # should always give norain if the threshold is set to 100%
-    assert check_norain(precip_arr, norain_thr=1.0)
+    assert check_norain(precip_arr, norain_thr=1.0, win_fun=None)
 
 
 # Finally, also test the compute_smooth_dilated mask functionality
diff --git a/pysteps/tests/test_nowcasts_anvil.py b/pysteps/tests/test_nowcasts_anvil.py
index 14a130fb1..48db86e60 100644
--- a/pysteps/tests/test_nowcasts_anvil.py
+++ b/pysteps/tests/test_nowcasts_anvil.py
@@ -1,3 +1,4 @@
+import numpy as np
 import pytest
 
 from pysteps import motion, nowcasts, verification
@@ -19,6 +20,28 @@
 ]
 
 
+def test_default_anvil_norain():
+    """Tests anvil nowcast with default params and all-zero inputs."""
+
+    # Define dummy nowcast input data
+    precip_input = np.zeros((4, 100, 100))
+
+    pytest.importorskip("cv2")
+    oflow_method = motion.get_method("LK")
+    retrieved_motion = oflow_method(precip_input)
+
+    nowcast_method = nowcasts.get_method("anvil")
+    precip_forecast = nowcast_method(
+        precip_input,
+        retrieved_motion,
+        timesteps=3,
+    )
+
+    assert precip_forecast.ndim == 3
+    assert precip_forecast.shape[0] == 3
+    assert precip_forecast.sum() == 0.0
+
+
 @pytest.mark.parametrize(anvil_arg_names, anvil_arg_values)
 def test_anvil_rainrate(
     n_cascade_levels,
diff --git a/pysteps/tests/test_nowcasts_linda.py b/pysteps/tests/test_nowcasts_linda.py
index 2d5f03b71..237dba4f0 100644
--- a/pysteps/tests/test_nowcasts_linda.py
+++ b/pysteps/tests/test_nowcasts_linda.py
@@ -26,6 +26,32 @@
 ]
 
 
+def test_default_linda_norain():
+    """Tests linda nowcast with default params and all-zero inputs."""
+
+    # Define dummy nowcast input data
+    precip_input = np.zeros((3, 100, 100))
+
+    pytest.importorskip("cv2")
+    oflow_method = motion.get_method("LK")
+    retrieved_motion = oflow_method(precip_input)
+
+    nowcast_method = nowcasts.get_method("linda")
+    precip_forecast = nowcast_method(
+        precip_input,
+        retrieved_motion,
+        n_ens_members=3,
+        timesteps=3,
+        kmperpixel=1,
+        timestep=5,
+    )
+
+    assert precip_forecast.ndim == 4
+    assert precip_forecast.shape[0] == 3
+    assert precip_forecast.shape[1] == 3
+    assert precip_forecast.sum() == 0.0
+
+
 @pytest.mark.parametrize(linda_arg_names, linda_arg_values)
 def test_linda(
     add_perturbations,
diff --git a/pysteps/tests/test_nowcasts_sprog.py b/pysteps/tests/test_nowcasts_sprog.py
index 1077c3edd..5872740e5 100644
--- a/pysteps/tests/test_nowcasts_sprog.py
+++ b/pysteps/tests/test_nowcasts_sprog.py
@@ -1,5 +1,6 @@
 # -*- coding: utf-8 -*-
 
+import numpy as np
 import pytest
 
 from pysteps import motion, nowcasts, verification
@@ -24,6 +25,29 @@
 ]
 
 
+def test_default_sprog_norain():
+    """Tests SPROG nowcast with default params and all-zero inputs."""
+
+    # Define dummy nowcast input data
+    precip_input = np.zeros((3, 100, 100))
+
+    pytest.importorskip("cv2")
+    oflow_method = motion.get_method("LK")
+    retrieved_motion = oflow_method(precip_input)
+
+    nowcast_method = nowcasts.get_method("sprog")
+    precip_forecast = nowcast_method(
+        precip_input,
+        retrieved_motion,
+        timesteps=3,
+        precip_thr=0.1,
+    )
+
+    assert precip_forecast.ndim == 3
+    assert precip_forecast.shape[0] == 3
+    assert precip_forecast.sum() == 0.0
+
+
 @pytest.mark.parametrize(sprog_arg_names, sprog_arg_values)
 def test_sprog(
     n_cascade_levels, ar_order, probmatching_method, domain, timesteps, min_csi
diff --git a/pysteps/tests/test_nowcasts_sseps.py b/pysteps/tests/test_nowcasts_sseps.py
index 4d89fd33a..b5ed73e6f 100644
--- a/pysteps/tests/test_nowcasts_sseps.py
+++ b/pysteps/tests/test_nowcasts_sseps.py
@@ -1,5 +1,6 @@
 # -*- coding: utf-8 -*-
 
+import numpy as np
 import pytest
 
 from pysteps import motion, nowcasts, verification
@@ -22,6 +23,37 @@
 ]
 
 
+def test_default_sseps_norain():
+    """Tests SSEPS nowcast with default params and all-zero inputs."""
+
+    # Define dummy nowcast input data
+    precip_input = np.zeros((3, 100, 100))
+    metadata = {
+        "accutime": 5,
+        "xpixelsize": 1000,
+        "threshold": 0.1,
+        "zerovalue": 0,
+    }
+
+    pytest.importorskip("cv2")
+    oflow_method = motion.get_method("LK")
+    retrieved_motion = oflow_method(precip_input)
+
+    nowcast_method = nowcasts.get_method("sseps")
+    precip_forecast = nowcast_method(
+        precip_input,
+        metadata,
+        retrieved_motion,
+        n_ens_members=3,
+        timesteps=3,
+    )
+
+    assert precip_forecast.ndim == 4
+    assert precip_forecast.shape[0] == 3
+    assert precip_forecast.shape[1] == 3
+    assert precip_forecast.sum() == 0.0
+
+
 @pytest.mark.parametrize(sseps_arg_names, sseps_arg_values)
 def test_sseps(
     n_ens_members,
diff --git a/pysteps/tests/test_nowcasts_steps.py b/pysteps/tests/test_nowcasts_steps.py
index 61af86ba5..7e558db45 100644
--- a/pysteps/tests/test_nowcasts_steps.py
+++ b/pysteps/tests/test_nowcasts_steps.py
@@ -7,7 +7,6 @@
 from pysteps import io, motion, nowcasts, verification
 from pysteps.tests.helpers import get_precipitation_fields
 
-
 steps_arg_names = (
     "n_ens_members",
     "n_cascade_levels",
@@ -22,7 +21,7 @@
 steps_arg_values = [
     (5, 6, 2, None, None, "spatial", 3, 1.30),
     (5, 6, 2, None, None, "spatial", [3], 1.30),
-    (5, 6, 2, "incremental", None, "spatial", 3, 7.31),
+    (5, 6, 2, "incremental", None, "spatial", 3, 7.32),
     (5, 6, 2, "sprog", None, "spatial", 3, 8.4),
     (5, 6, 2, "obs", None, "spatial", 3, 8.37),
     (5, 6, 2, None, "cdf", "spatial", 3, 0.60),
@@ -31,6 +30,33 @@
 ]
 
 
+def test_default_steps_norain():
+    """Tests STEPS nowcast with default params and all-zero inputs."""
+
+    # Define dummy nowcast input data
+    precip_input = np.zeros((3, 100, 100))
+
+    pytest.importorskip("cv2")
+    oflow_method = motion.get_method("LK")
+    retrieved_motion = oflow_method(precip_input)
+
+    nowcast_method = nowcasts.get_method("steps")
+    precip_forecast = nowcast_method(
+        precip_input,
+        retrieved_motion,
+        n_ens_members=3,
+        timesteps=3,
+        precip_thr=0.1,
+        kmperpixel=1,
+        timestep=5,
+    )
+
+    assert precip_forecast.ndim == 4
+    assert precip_forecast.shape[0] == 3
+    assert precip_forecast.shape[1] == 3
+    assert precip_forecast.sum() == 0.0
+
+
 @pytest.mark.parametrize(steps_arg_names, steps_arg_values)
 def test_steps_skill(
     n_ens_members,
diff --git a/pysteps/utils/check_norain.py b/pysteps/utils/check_norain.py
new file mode 100644
index 000000000..f4e3c7058
--- /dev/null
+++ b/pysteps/utils/check_norain.py
@@ -0,0 +1,51 @@
+import numpy as np
+
+from pysteps import utils
+
+
+def check_norain(precip_arr, precip_thr=None, norain_thr=0.0, win_fun=None):
+    """
+
+    Parameters
+    ----------
+    precip_arr:  array-like
+      An at least 2 dimensional array containing the input precipitation field
+    precip_thr: float, optional
+      Specifies the threshold value for minimum observable precipitation intensity. If None, the
+      minimum value over the domain is taken.
+    norain_thr: float, optional
+      Specifies the threshold value for the fraction of rainy pixels in precip_arr below which we consider there to be
+      no rain. Standard set to 0.0
+    win_fun: {'hann', 'tukey', None}
+      Optional tapering function to be applied to the input field, generated with
+      :py:func:`pysteps.utils.tapering.compute_window_function`
+      (default None).
+      This parameter needs to match the window function you use in later noise generation,
+      or else this method will say that there is rain, while after the tapering function is
+      applied there is no rain left, so you will run into a ValueError.
+    Returns
+    -------
+    norain: bool
+      Returns whether the fraction of rainy pixels is below the norain_thr threshold.
+
+    """
+
+    if win_fun is not None:
+        tapering = utils.tapering.compute_window_function(
+            precip_arr.shape[-2], precip_arr.shape[-1], win_fun
+        )
+    else:
+        tapering = np.ones((precip_arr.shape[-2], precip_arr.shape[-1]))
+
+    tapering_mask = tapering == 0.0
+    masked_precip = precip_arr.copy()
+    masked_precip[..., tapering_mask] = np.nanmin(precip_arr)
+
+    if precip_thr is None:
+        precip_thr = np.nanmin(masked_precip)
+    rain_pixels = masked_precip[masked_precip > precip_thr]
+    norain = rain_pixels.size / masked_precip.size <= norain_thr
+    print(
+        f"Rain fraction is: {str(rain_pixels.size / masked_precip.size)}, while minimum fraction is {str(norain_thr)}"
+    )
+    return norain
diff --git a/pysteps/utils/tapering.py b/pysteps/utils/tapering.py
index aeb4100eb..b4073ec9e 100644
--- a/pysteps/utils/tapering.py
+++ b/pysteps/utils/tapering.py
@@ -82,7 +82,7 @@ def compute_window_function(m, n, func, **kwargs):
         Array of shape (m, n) containing the tapering weights.
     """
     X, Y = np.meshgrid(np.arange(n), np.arange(m))
-    R = np.sqrt((X - int(n / 2)) ** 2 + (Y - int(m / 2)) ** 2)
+    R = np.sqrt(((X / n) - 0.5) ** 2 + ((Y / m) - 0.5) ** 2)
 
     if func == "hann":
         return _hann(R)
@@ -108,26 +108,24 @@ def _compute_mask_distances(mask):
 
 def _hann(R):
     W = np.ones_like(R)
-    N = min(R.shape[0], R.shape[1])
-    mask = R > int(N / 2)
+    mask = R > 0.5
 
     W[mask] = 0.0
-    W[~mask] = 0.5 * (1.0 - np.cos(2.0 * np.pi * (R[~mask] + int(N / 2)) / N))
+    W[~mask] = 0.5 * (1.0 - np.cos(2.0 * np.pi * (R[~mask] + 0.5)))
 
     return W
 
 
 def _tukey(R, alpha):
     W = np.ones_like(R)
-    N = min(R.shape[0], R.shape[1])
 
-    mask1 = R < int(N / 2)
-    mask2 = R > int(N / 2) * (1.0 - alpha)
+    mask1 = R < 0.5
+    mask2 = R > 0.5 * (1.0 - alpha)
     mask = np.logical_and(mask1, mask2)
     W[mask] = 0.5 * (
-        1.0 + np.cos(np.pi * (R[mask] / (alpha * 0.5 * N) - 1.0 / alpha + 1.0))
+        1.0 + np.cos(np.pi * (R[mask] / (alpha * 0.5) - 1.0 / alpha + 1.0))
     )
-    mask = R >= int(N / 2)
+    mask = R >= 0.5
     W[mask] = 0.0
 
     return W