Repository: pySTEPS/pysteps

Issue:
Initializing noise field with default arguments returns array of zeros. Is this expected behavior?
I am running a probabilistic nowcast for a rectangular precipitation field (199x113) using the steps function with the following settings: 
```
domain = 'spatial'
extrap_method = 'semilagrangian'
decomp_method = 'fft'
fft_method = 'numpy'
noise_method = 'nonparametric'
noise_stddev_adj = 'auto'
```
 and run into the following error before the noise adjustment coefficients are computed: 
```ValueError: field contains non-finite values```. 

This points me to the following statement in the cascade [decomposition](https://github.com/pySTEPS/pysteps/blob/e3325854dae006b1e78eb43a6d2203f2b0b71560/pysteps/cascade/decomposition.py#L193C1-L194C61):
```
if np.any(~np.isfinite(field)):
    raise ValueError("field contains non-finite values") 
```

which in turn points me to ``` decomp_N = decomp_method(N, F, mask=MASK_)``` in the following function: https://github.com/pySTEPS/pysteps/blob/e3325854dae006b1e78eb43a6d2203f2b0b71560/pysteps/noise/utils.py#L24 

This stems from the fact that the noise field F returned in the function:
https://github.com/pySTEPS/pysteps/blob/e3325854dae006b1e78eb43a6d2203f2b0b71560/pysteps/noise/fftgenerators.py#L214

is an array of zeroes, and thus N in the following lines in [noise generation](https://github.com/pySTEPS/pysteps/blob/e3325854dae006b1e78eb43a6d2203f2b0b71560/pysteps/noise/fftgenerators.py#L400C1-L418C75) returns a ```nan```. 

It appears the issue lies in the following lines in the function https://github.com/pySTEPS/pysteps/blob/e3325854dae006b1e78eb43a6d2203f2b0b71560/pysteps/noise/fftgenerators.py#L214
```
F = np.zeros(fft_shape, dtype=complex)
 for i in range(nr_fields):
     if use_full_fft:
         F += fft.fft2(field[i, :, :] * tapering)
     else:
         F += fft.rfft2(field[i, :, :] * tapering)
 F /= nr_fields
```

Using the default arguments, with window_function="tukey" and a field size of 199x113, the tapering returns an oval shape (black for values >0):  
![download](https://github.com/user-attachments/assets/6c749cb8-a2dd-4a22-a265-9bb56571709a)

This essentially means that if we have a precipitation field, with precipitation only at the edges, and use the default steps settings, this would return a noise field of zeros and eventually lead to an error in the cascade decomposition. 

Though the function appears to numerically behave as expected, is this an error or is this desired behavior? Does this mean it's not possible to compute a nowcast if we have precipitation only around the edges of our domain?

Discussion:
@bwalraven: Adding a small piece of code as example: 
```
# Initialize zero array
rainy_edges = np.zeros((3, 199, 113))

# Disturb the top 20 and bottom 20 rows for each field with values between 1-10
for i in np.arange(0, 3):
    rainy_edges[i, :20, :] = np.random.randint(1, 10, (20, 113))
    rainy_edges[i, -20:, :] = np.random.randint(1, 10, (20, 113))


# Set default arguments
field = rainy_edges
use_full_fft = 'False' # default
win_fun = 'tukey' # default

# Copied code from the non-parametric noise initializer
# dims
if len(field.shape) == 2:
    field = field[None, :, :]
nr_fields = field.shape[0]
field_shape = field.shape[1:]
if use_full_fft:
    fft_shape = (field.shape[1], field.shape[2])
else:
    fft_shape = (field.shape[1], int(field.shape[2] / 2) + 1)

# make sure non-rainy pixels are set to zero
field -= field.min(axis=(1, 2))[:, None, None]

if win_fun is not None:
    tapering = utils.tapering.compute_window_function(
        field_shape[0], field_shape[1], win_fun
    )
else:
    tapering = np.ones(field_shape)

F = np.zeros(fft_shape, dtype=complex)
for i in range(nr_fields):
    if use_full_fft:
        F += fft.fft2(field[i, :, :] * tapering)
    else:
        F += fft.rfft2(field[i, :, :] * tapering)
F /= nr_fields

# returns a field with only zeros
print(F.real.sum())
```

```
# plotting the tapering
plt.figure(figsize=(6, 10))
plt.title('Tapering')
plt.pcolormesh(tapering>0, cmap='Greys')
````@dnerini: Hi @bwalraven thanks for the detailed issue report. 

I'll try to have a close look, but in essence the tapering is probably behaving as expected as we exclude rainy pixels at the edges when computing the power spectra to avoid artefacts. 

However, the plot of the tapering function does look a bit weird, and, more importantly, it should handle zero fields more graciously instead of raising an error. @dnerini: I also wonder if this is not an edge case that was not covered when addressing #309 @ladc @bwalraven: Hi @dnerini thanks for looking into this. Curious to what you can find out. 

And understandable to exclude some pixels at the edges, though what strikes me in the example I showed is the difference in excluded pixels in x and y direction (approx. 1 pixel vs. almost 50 pixels). Note that the image I included is a little bit distorted as there are more y-pixels (199) than x-pixels (113).  

@mats-knmi: We get the same errors semi-regularly, when there is only a small amount of precipitation in the model. There are 2 issues I see here:
1. The tapering method
For us the tapering does seem to look alright (for a domain of 700x765 pixels). It looks to me like this "tukey" tapering method always procuces a perfect circle that fits the domain. In the case of @bwalraven (for a very rectangular domain) this is not a good tapering method, and I think this should maybe not be the default. Maybe we can stretch this tapering to match the input field shape (and produce an oval tapering).
![Image](https://github.com/user-attachments/assets/dc9bcde6-e0a7-4204-bc15-3dc20ad6ff60)
2. The error handling
In the solution of issue #309 some checks were implemented for no rain in NWP and/or QPE. However, these checks are done before the tapering is applied. This results in the issue described here. In my opinion we should check if applying the tapering to `self.__precip` results in zero precip and then handle that the same as `self.__params.zero_precip_radar`. Next we should check if applying the tapering to `self.__state.precip_noise_input` results in zero precip and then handle that the same as `self.__params.zero_precip_model_fields` (i.e. returning `self.__zero_precipitation_forecast()`). I will create a PR with this proposed solution. Let me know what you think @dnerini. (but also @ladc and @RubenImhoff, since you worked on #309).@RubenImhoff: Thanks, @mats-knmi! I think combining both points that you mention is a good way to find a solution for this issue. :) Happy to give the proposed solution a review.@mats-knmi: I created a PR: https://github.com/pySTEPS/pysteps/pull/460@SKB-7: I'm also facing the same issue as mentioned in #451, " is an array of zeroes, and thus N in the following lines in [noise generation] returns a nan."

I have seen PR: #460 created by @mats-knmi. It is not completely merged yet, right? Does something need to be fixed? I'm happy to get involved.

Many Thanks.@mats-knmi: Hi @SKB-7, thanks for your enthusiasm. Technically #460 is ready to be merged. I am just waiting for @dnerini to answer, since he first picked up this issue.

It would be helpful for us to know if that PR solves your issue though. Can you try it out and see if it works? You could try to do the following to test that PR out.
```
pip install git+https://github.com/pySTEPS/pysteps.git@fix-noise-tapering-and-error-handling
```@SKB-7: Hi @mats-knmi, After installing the branch pysteps.git@fix-noise-tapering-and-error-handling, I have tested my code, but still, noise generation using initialize_nonparam_2d_ssft_filter and generate_noise_2d_ssft_filter returns a nan values.@mats-knmi: @SKB-7 Interesting, which part of pysteps do you run? My PR should fix this issue for the steps blending. Do you maybe have an example of how you run pysteps, so that I can reproduce this?@SKB-7: Below is the function I'm using to generate noise for n ensembles.

def getCorrNoise(n, h, dt, obsFile, U_windFile, V_windFile, seednum):

    ds1 = xr.open_dataset(U_windFile)
    ds2 = xr.open_dataset(V_windFile)
    
    u_values = ds1['u'][0:(0+n),:,:] # eastward wind [m/s]
    v_values = ds2['v'][0:(0+n),:,:] # northward wind [m/s]
    
    # ----  open precip data  ----
    ds = xr.open_dataset(obsFile)
    
    har = ds['prcp'][0:(0+n),:,:].data
    har[har<0.1] = 0 
    ds=None
    xsize=har.shape[2]
    ysize=har.shape[1]
    
    yind=np.repeat(np.arange(0,ysize),xsize).reshape(ysize,xsize).astype('int16')
    xind=np.tile(np.arange(0,xsize),ysize).reshape(ysize,xsize).astype('int16')
 
    dx = dy = 10000
    
    # create array to hold simulated noise
    s = np.empty((n,ysize,xsize),dtype=np.float32)

    # --- SIMULATE FIRST FIELD ---
    if len(har[0,:,:][har[0,:,:]>=0.1]) > 1133:          
        # print("entered 1st if")
        # only use pysteps smoothing filter if there is actually rain in the study area
        Fnp = initialize_nonparam_2d_ssft_filter(har[0,:,:])#, win_size=(64,64)) 
        s[0,:,:] = generate_noise_2d_ssft_filter(Fnp, seed=seednum)
        firstRain=True
        
    else:
        # print("entered 1st else")
        # otherwise just use smoothed white noise
        s[0,:,:] = np.random.normal(size=(ysize,xsize)) 
        #**seed kwargs are not there in normal fnc so commented,seed=seednum)
        firstRain=False
    
    last_rn = s[0,:,:] # record this instance of random noise as most recent instance of random noise
    
    
    # set initial value of alpha parameter - determines degree to which noise is incorporated into field at each timestep
    alpha = 0.95
    
    # --- ADVECT FIELD FORWARD IN TIME ---
    for hr in range(1,n):
        # print(f"hr: {hr}")
        
        # get wind values from PREVIOUS DAY and calculate dx, dy that correspond to this time step
        u = u_values[hr-1,:,:] # eastward wind [m/s]
        deltax = u*60.*60.*24. # distance wind travels each time step in positive x dir [m]
        dix = deltax//dx # number of grid steps to shift in semi-lagrangian scheme in units of [0.1 deg pixels]
        
        v = v_values[hr-1,:,:] # northward wind [m/s]
        deltay = -v*60.*60.*24. # distance wind travels each time step in positive y dir [m]
        diy = deltay//dy # number of grid steps to shift in semi-lagrangian scheme [0.1 deg pixels]
        
        harh = har[hr,:,:] 
        
        # use pySTEPS to generate noise with power spectrum of Precip field
        # only use pysteps smoothing filter if there is actually rain in the study area
        
        if len(harh[harh>=0.1]) > 1133: # atleast 5% of study area must be rainy
            # print("entered 2nd if")
            
            Fnp = initialize_nonparam_2d_ssft_filter(harh)#, win_size=(64,64))
            rn = generate_noise_2d_ssft_filter(Fnp, seed=seednum+hr)
            # print(f'seednum {seednum+hr}')
            # get average temporal correlation of field in previous timestep
            alpha = getAlpha(harh,har[hr-1,:,:])
           
            firstRain=True            
        else:
            # print("entered 2nd else")
            
            # if field does not have at least 5% rainy pixels, generate random field 
            # using spatial correlation structure of last instance of rainfall 
            
            if firstRain==True:
                print("entered 3rd if")                
                rn = generate_noise_2d_ssft_filter(Fnp, seed=seednum+hr)
            
            else:
                # print("entered 3rd else")
                # if first instance of a "rainy field" hasn't occurred yet in study period, use white noise
                rn = np.random.normal(size=(ysize,xsize))
                # **seed kwargs are not there in normal fnc so commented,seed=seednum+hr) # white noise
                # print(f'random noise {hr}')
                        
        
        if np.isnan(rn[0,0])==True: # if pysteps smoothing process doesn't work for some reason, use last realization of smoothed noise
            print("NaN in noise field at timestep %d, using last realization of noise field"%hr)
            rn = last_rn
            # print(f'error ssft_filter {hr}')
        else:
            last_rn = rn
        
    
        # -- advect noise values from previous time step & perturb with noise field rn --
        ybefore=np.array((yind-diy)%ysize,dtype='int16')
        xbefore=np.array((xind-dix)%xsize,dtype='int16')
        s[hr,:] = (alpha*s[hr-1,ybefore,xbefore] + np.sqrt(1.-alpha**2)*rn).astype('float32')
        # print(f"hr: {hr} ,Size of s: {s.shape}")
    
    # return noise field
    return(s)
@mats-knmi: @SKB-7 I think I understand the issue now. I have fixed this for the [steps blending](https://github.com/pySTEPS/pysteps/blob/master/pysteps/blending/steps.py), not for any of the nowcasts (and also not if you call these noise generation functions directly, like you are). Re-reading the issue I see that @bwalraven probably also had the issue with the nowcasts and not the blending.

In these cases there is currently no check for no rain built into the methods themselves, so you have to write them yourself before calling the method. I see that you have already implemented this in your code `if len(harh[harh>=0.1]) > 1133`, but this does not cover the case where there is only precipitation at the edges.

Basically the way I solved this for the blending is by improving on that check. I suppose a solution would be to add this improved `norain` check as a util function to the pysteps code, so that you can use this to check if you can call pysteps noise generation, so that you can decide on a fallback.

Also, just out of curiosity, what are you using pysteps for? If you are generating an ensemble nowcast, why are you not using the [steps nowcast](https://github.com/pySTEPS/pysteps/blob/master/pysteps/nowcasts/steps.py) (or any of the other [nowcast modules](https://github.com/pySTEPS/pysteps/tree/master/pysteps/nowcasts))? And if you are creating your own improved nowcast method, would you maybe consider contributing your alternative method to the pysteps project?@mats-knmi: @SKB-7 I created a new [check_norain](https://github.com/pySTEPS/pysteps/blob/fix-noise-tapering-and-error-handling/pysteps/utils/check_norain.py#L6) method, you should be able to use this in stead of your if statement and then you should not be getting any errors anymore. Please let me know if this works for you!@bwalraven: > Re-reading the issue I see that [@bwalraven](https://github.com/bwalraven) probably also had the issue with the nowcasts and not the blending.

You are right @mats-knmi. I only use it for generating ensemble nowcasts and don't do any blending. 
I use a custom `check_norain` function (much like the one you just added to the `pysteps.utils`, thanks!), and in my current workflow have also added a `check_zero_array_in_noise_generation` function before starting a nowcast. If any of the two returns `True` the nowcast is not run but simply returns a field of zeroes. 

```
def check_zero_array_in_noise_generation(
    field: npt.NDArray,
    win_fun: str = 'tukey',
    use_full_fft: bool = False,
    fft: str = 'numpy'
) -> bool:
    """
    Function to check if the noise generator field contains only zeroes after multiplication with the tapering.  
    This yields an error a.o. in the cascade decomposition. 
    
    Parameters
    ----------
    field:  array-like
      Array containing the input precipitation field
    win_fun: str, optional
      Tapering window function. Default is 'tukey'. 
    use_full_fft: bool, optional
      Whether to use the full fast fourier transform. Default is False.
    fft: str, optional
       Fast Fourier Transform method. Default 'numpy'.  
      
    Returns
    -------
    zeros_array: bool
      Returns True if the generated noise field contains only zero values, otherwise False.
    """
    # replace non-finite values with the minimum value
    field = field.copy()    
    for i in range(field.shape[0]):
        field[i, ~np.isfinite(field[i, :])] = np.nanmin(field[i, :])


    fft_shape = field.shape if len(field.shape) == 2 else field.shape[1:]
    fft = pysteps.utils.get_method(fft, shape=fft_shape)

    nr_fields = field.shape[0]
    field_shape = field.shape[1:]
    if use_full_fft:
        fft_shape = (field.shape[1], field.shape[2])
    else:
        fft_shape = (field.shape[1], int(field.shape[2] / 2) + 1)
    
    # make sure non-rainy pixels are set to zero
    field -= field.min(axis=(1, 2))[:, None, None]
    
    if win_fun is not None:
        tapering = pysteps.utils.tapering.compute_window_function(
            field_shape[0], field_shape[1], win_fun
        )
    else:
        tapering = np.ones(field_shape)
    
    F = np.zeros(fft_shape, dtype=complex)
    for i in range(nr_fields):
        if use_full_fft:
            F += fft.fft2(field[i, :, :] * tapering)
        else:
            F += fft.rfft2(field[i, :, :] * tapering)
    F /= nr_fields

    zeros_array = np.all(F == 0)
    if zeros_array:
        print("Generated noise field contains only zeros")

    return zeros_array
```@dnerini: thanks @bwalraven for sharing your thoughts. Why do you return fields of zeros in this case? Would it make more sense to abort computation of the nowcast instead?@dnerini: Looks to me that #460 should include something that addresses this issue for the nowcast only method, would you agree @mats-knmi ? We could take a simpler approach in this case, as raising a more informative error message if all fields passed to initialization of the noise generator are zeros@bwalraven: > thanks [@bwalraven](https://github.com/bwalraven) for sharing your thoughts. Why do you return fields of zeros in this case? Would it make more sense to abort computation of the nowcast instead?

Yes, essentially that is what I meant. Since I am running the nowcasts for many timesteps, if for a timestep the `check_norain` or `check_zero_array_in_noise_generation` is not passed, it returns and empty nowcasts, and moves on to the next timestep. But I still save the empty nowcast to file for a smooth workflow down the road. @mats-knmi: @dnerini Yes, I initially missed that side of this issue. But I think that the current solution actually handles this very well, because users can use the `check_norain` method before starting any nowcast to make sure they don't run into this error. I will also update the error message to tell the users that this is what they need to do.

Finally I also believe it could be beneficial to call the `check_norain` method from the nowcast methods themselves (as we already do for the blending) and then return a nowcast with only zeros if the check fails. This would completely remove the issue from the user side, except for a case like @SKB-7 where you call the fftgenerators directly, but then the updated error message would still help.@SKB-7: > [@SKB-7](https://github.com/SKB-7) I created a new [check_norain](https://github.com/pySTEPS/pysteps/blob/fix-noise-tapering-and-error-handling/pysteps/utils/check_norain.py#L6) method, you should be able to use this in stead of your if statement and then you should not be getting any errors anymore. Please let me know if this works for you!

Hi @mats-knmi, sorry for the late reply. I've tested with my code and data. If rain_pixels.size is 0 in the check_norain.py file then using initialize_nonparam_2d_ssft_filter rises a "ValueError: zero-size array to reduction operation minimum which has no identity".@mats-knmi: @SKB-7 If `rain_pixels.size` is `0`, then `check_norain` should return `True` right? So then you should not call `initialize_nonparam_2d_ssft_filter`, since that method only works when there is rain.@SKB-7: @mats-knmi Yeah, right. Sorry, Just got confused. I was getting the following error:


File ~/envs/pysteps_env/lib/python3.10/site-packages/pysteps/noise/fftgenerators.py:509, in initialize_nonparam_2d_ssft_filter(field, **kwargs)
    507 # remove rain/no-rain discontinuity
    508 if rm_rdisc:
--> 509     field[field > field.min()] -= field[field > field.min()].min() - field.min()
    511 # dims
    512 if len(field.shape) == 2:

File ~/envs/pysteps_env/lib/python3.10/site-packages/numpy/_core/_methods.py:48, in _amin(a, axis, out, keepdims, initial, where)
     46 def _amin(a, axis=None, out=None, keepdims=False,
     47           initial=_NoValue, where=True):
---> 48     return umr_minimum(a, axis, None, out, keepdims, initial, where)

ValueError: zero-size array to reduction operation minimum which has no identity@ladc: Hi everyone, I'm a bit late to the conversation but indeed the no-rain-check was only implemented for the blending and not generally for all nowcast methods. Possibly this can be checked as a step in all nowcast methods, as suggested by @mats-knmi? The `win_fun` can then be filled in with a sane value depending on the method used instead of leaving this to the user (susceptible to more errors)?

Solution patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index cb599f2..cb9b634 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,6 +1,6 @@
 repos:
 -   repo: https://github.com/psf/black
-    rev: 24.10.0
+    rev: 25.1.0
     hooks:
     - id: black
       language_version: python3
diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py
index 740fce8..c722e43 100644
--- a/pysteps/blending/steps.py
+++ b/pysteps/blending/steps.py
@@ -56,6 +56,7 @@ from pysteps import blending, cascade, extrapolation, noise, utils
 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 @@ class StepsBlendingNowcaster:
 
         # 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 @@ class StepsBlendingNowcaster:
             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 @@ class StepsBlendingNowcaster:
             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 @@ class StepsBlendingNowcaster:
                 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 @@ class StepsBlendingNowcaster:
             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 @@ class StepsBlendingNowcaster:
             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 486d7af..aaed2cf 100644
--- a/pysteps/blending/utils.py
+++ b/pysteps/blending/utils.py
@@ -15,11 +15,11 @@ Module with common utilities used by the blending methods.
     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 import get_method as cascade_get_method
 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/datasets.py b/pysteps/datasets.py
index b90ea0b..91abc62 100644
--- a/pysteps/datasets.py
+++ b/pysteps/datasets.py
@@ -21,7 +21,6 @@ import shutil
 import sys
 import time
 from datetime import datetime, timedelta
-from distutils.dir_util import copy_tree
 from logging.handlers import RotatingFileHandler
 from tempfile import NamedTemporaryFile, TemporaryDirectory
 from urllib import request
@@ -329,7 +328,9 @@ def download_pysteps_data(dir_path, force=True):
 
         zip_obj.extractall(tmp_dir.name)
 
-        copy_tree(os.path.join(tmp_dir.name, common_path), dir_path)
+        shutil.copytree(
+            os.path.join(tmp_dir.name, common_path), dir_path, dirs_exist_ok=True
+        )
 
 
 def create_default_pystepsrc(
diff --git a/pysteps/noise/fftgenerators.py b/pysteps/noise/fftgenerators.py
index 7ffabdc..d3414d0 100644
--- a/pysteps/noise/fftgenerators.py
+++ b/pysteps/noise/fftgenerators.py
@@ -46,6 +46,7 @@ field of correlated noise cN of shape (m, n).
 
 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 5bd35e1..cee9902 100644
--- a/pysteps/nowcasts/linda.py
+++ b/pysteps/nowcasts/linda.py
@@ -40,6 +40,8 @@ the short-space Fourier transform (SSFT) methodology developed in
 import time
 import warnings
 
+from pysteps.utils.check_norain import check_norain
+
 try:
     import dask
 
@@ -47,14 +49,14 @@ try:
 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 3d9e78f..3742556 100644
--- a/pysteps/nowcasts/sprog.py
+++ b/pysteps/nowcasts/sprog.py
@@ -10,16 +10,16 @@ Implementation of the S-PROG method described in :cite:`Seed2003`
     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 94cd345..f2891b5 100644
--- a/pysteps/nowcasts/sseps.py
+++ b/pysteps/nowcasts/sseps.py
@@ -18,17 +18,16 @@ size.
     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 818123d..dc77c7e 100644
--- a/pysteps/nowcasts/steps.py
+++ b/pysteps/nowcasts/steps.py
@@ -11,22 +11,24 @@ Implementation of the STEPS stochastic nowcasting method as described in
     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 @@ class StepsNowcaster:
         # 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 @@ class StepsNowcaster:
             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 fd111e2..8ddd3da 100644
--- a/pysteps/nowcasts/utils.py
+++ b/pysteps/nowcasts/utils.py
@@ -17,6 +17,7 @@ Module with common utilities used by nowcasts methods.
 """
 
 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/utils/check_norain.py b/pysteps/utils/check_norain.py
new file mode 100644
index 0000000..f4e3c70
--- /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 aeb4100..b4073ec 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


Test patch:
diff --git a/pysteps/tests/test_blending_utils.py b/pysteps/tests/test_blending_utils.py
index 22d0a00..401b6f1 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 14a130f..48db86e 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 @@ anvil_arg_values = [
 ]
 
 
+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 2d5f03b..237dba4 100644
--- a/pysteps/tests/test_nowcasts_linda.py
+++ b/pysteps/tests/test_nowcasts_linda.py
@@ -26,6 +26,32 @@ linda_arg_values = [
 ]
 
 
+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 1077c3e..5872740 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 @@ sprog_arg_values = [
 ]
 
 
+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 4d89fd3..b5ed73e 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 @@ sseps_arg_values = [
 ]
 
 
+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 61af86b..7e558db 100644
--- a/pysteps/tests/test_nowcasts_steps.py
+++ b/pysteps/tests/test_nowcasts_steps.py
@@ -7,7 +7,6 @@ import pytest
 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_names = (
 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 @@ steps_arg_values = [
 ]
 
 
+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,