fix noise tapering and error handling

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
@@ -777,7 +778,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 +1094,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):
@@ -1177,10 +1180,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[

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 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 @@
       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 in stead"
     )
-    return norain
+    return new_check_norain(precip_arr, precip_thr, norain_thr, None)
 
 
 def compute_smooth_dilated_mask(

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

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),

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="tukey"):
+    """
+
+    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 'tukey').
+      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

pysteps/utils/tapering.py

@@ -82,7 +82,7 @@
         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 _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