Fix warnings shown during tests (#152)

* pysteps.motion.vet: Fix bug in `vet` keywords handling. Previously, the argument passed to the minimizer was ignored.
* pysteps.io.importers: Add a function to safely get the rain/no rain threshold that handles NaN correctly.
* Fix collections imports. Remove DeprecationWarning for using the ABCs from 'collections' instead  of from 'collections.abc'.
* pysteps.timeseries.autoregression: Make some docstrings raw literals, avoiding the interpretation of math equations as invalid escape sequences.

Author: aperezhortal

pysteps/io/importers.py

@@ -201,6 +201,31 @@ def _get_grib_projection(grib_msg):
     return projparams
 
 
+def _get_threshold_value(precip):
+    """
+    Get the the rain/no rain threshold with the same unit, transformation and
+    accutime of the data.
+    If all the values are NaNs, the returned value is `np.nan`.
+    Otherwise, np.min(precip[precip > precip.min()]) is returned.
+
+    Returns
+    -------
+
+    threshold : float
+    """
+    valid_mask = np.isfinite(precip)
+    if valid_mask.any():
+        _precip = precip[valid_mask]
+        min_precip = _precip.min()
+        above_min_mask = _precip > min_precip
+        if above_min_mask.any():
+            return np.min(_precip[above_min_mask])
+        else:
+            return min_precip
+    else:
+        return np.nan
+      
+
 @postprocess_import(dtype='float32')
 def import_mrms(filename, extent=None, window_size=4, **kwargs):
     """
@@ -371,6 +396,7 @@ def import_mrms(filename, extent=None, window_size=4, **kwargs):
         zerovalue=0,
         projection=proj_params,
         yorigin="upper",
+        threshold=_get_threshold_value(precip),
         x1=x1,
         x2=x2,
         y1=y1,
@@ -414,10 +440,7 @@ def import_bom_rf3(filename, **kwargs):
 
     metadata["transform"] = None
     metadata["zerovalue"] = np.nanmin(precip)
-    if np.any(np.isfinite(precip)):
-        metadata["threshold"] = np.nanmin(precip[precip > np.nanmin(precip)])
-    else:
-        metadata["threshold"] = np.nan
+    metadata["threshold"] = _get_threshold_value(precip)
 
     return precip, None, metadata
 
@@ -584,9 +607,8 @@ def import_fmi_geotiff(filename, **kwargs):
     metadata["unit"] = rb.GetUnitType()
     metadata["transform"] = None
     metadata["accutime"] = 5.0
-    precip_min = np.nanmin(precip)
-    metadata["threshold"] = np.nanmin(precip[precip > precip_min])
-    metadata["zerovalue"] = precip_min
+    metadata["threshold"] = _get_threshold_value(precip)
+    metadata["zerovalue"] = np.nanmin(precip)
 
     return precip, None, metadata
 
@@ -640,10 +662,7 @@ def import_fmi_pgm(filename, gzipped=False, **kwargs):
     metadata["unit"] = "dBZ"
     metadata["transform"] = "dB"
     metadata["zerovalue"] = np.nanmin(precip)
-    if np.any(np.isfinite(precip)):
-        metadata["threshold"] = np.nanmin(precip[precip > np.nanmin(precip)])
-    else:
-        metadata["threshold"] = np.nan
+    metadata["threshold"] = _get_threshold_value(precip)
     metadata["zr_a"] = 223.0
     metadata["zr_b"] = 1.53
 
@@ -861,7 +880,7 @@ def import_knmi_hdf5(filename, qty="ACRR", accutime=5.0, pixelsize=1.0, **kwargs
     metadata["unit"] = unit
     metadata["transform"] = transform
     metadata["zerovalue"] = 0.0
-    metadata["threshold"] = np.nanmin(precip[precip > np.nanmin(precip)])
+    metadata["threshold"] = _get_threshold_value(precip)
     metadata["zr_a"] = 200.0
     metadata["zr_b"] = 1.6
 
@@ -985,10 +1004,7 @@ def import_mch_gif(filename, product, unit, accutime, **kwargs):
     metadata["unit"] = unit
     metadata["transform"] = None
     metadata["zerovalue"] = np.nanmin(precip)
-    if np.any(precip > np.nanmin(precip)):
-        metadata["threshold"] = np.nanmin(precip[precip > np.nanmin(precip)])
-    else:
-        metadata["threshold"] = np.nan
+    metadata["threshold"] = _get_threshold_value(precip)
     metadata["institution"] = "MeteoSwiss"
     metadata["product"] = product
     metadata["zr_a"] = 316.0
@@ -1201,10 +1217,7 @@ def import_mch_metranet(filename, product, unit, accutime):
     metadata["unit"] = unit
     metadata["transform"] = None
     metadata["zerovalue"] = np.nanmin(precip)
-    if np.isnan(metadata["zerovalue"]):
-        metadata["threshold"] = np.nan
-    else:
-        metadata["threshold"] = np.nanmin(precip[precip > metadata["zerovalue"]])
+    metadata["threshold"] = _get_threshold_value(precip)
     metadata["zr_a"] = 316.0
     metadata["zr_b"] = 1.5
 
@@ -1389,11 +1402,6 @@ def import_opera_hdf5(filename, qty="RATE", **kwargs):
         unit = "mm/h"
         transform = None
 
-    if np.any(np.isfinite(precip)):
-        thr = np.nanmin(precip[precip > np.nanmin(precip)])
-    else:
-        thr = np.nan
-
     metadata = {
         "projection": proj4str,
         "ll_lon": ll_lon,
@@ -1412,7 +1420,7 @@ def import_opera_hdf5(filename, qty="RATE", **kwargs):
         "unit": unit,
         "transform": transform,
         "zerovalue": np.nanmin(precip),
-        "threshold": thr,
+        "threshold": _get_threshold_value(precip),
     }
 
     f.close()
@@ -1493,10 +1501,7 @@ def import_saf_crri(filename, extent=None, **kwargs):
 
     metadata["transform"] = None
     metadata["zerovalue"] = np.nanmin(precip)
-    if np.any(np.isfinite(precip)):
-        metadata["threshold"] = np.nanmin(precip[precip > np.nanmin(precip)])
-    else:
-        metadata["threshold"] = np.nan
+    metadata["threshold"] = _get_threshold_value(precip)
 
     return precip, quality, metadata
 

pysteps/motion/vet.py

@@ -481,6 +481,7 @@ def debug_print(*args, **kwargs):
     options.setdefault('gtol', 0.1)
     options.setdefault('maxiter', 100)
     options.setdefault('disp', False)
+    optimization_method = options.pop("method", "CG")
 
     # Set to None to suppress pylint warning.
     pad_i = None
@@ -613,7 +614,7 @@ def debug_print(*args, **kwargs):
                                 (sectors_in_i, sectors_in_j),
                                 _mask,
                                 smooth_gain),
-                          method='CG',
+                          method=optimization_method,
                           options=options)
 
         first_guess = result.x.reshape(*first_guess.shape)

pysteps/tests/test_motion.py

@@ -209,7 +209,7 @@ def test_optflow_method_convergence(input_precip, optflow_method_name,
         # To increase the stability of the tests to we increase this value to
         # maxiter=150.
         retrieved_motion = oflow_method(precip_obs, verbose=False,
-                                       options=dict(maxiter=150, method='BFGS'))
+                                       options=dict(maxiter=150))
     elif optflow_method_name == 'proesmans':
         retrieved_motion = oflow_method(precip_obs)
     else:

pysteps/timeseries/autoregression.py

@@ -116,7 +116,7 @@ def ar_acf(gamma, n=None):
 
 def estimate_ar_params_ols(x, p, d=0, check_stationarity=True,
                            include_constant_term=False, h=0, lam=0.0):
-    """Estimate the parameters of an autoregressive AR(p) model
+    r"""Estimate the parameters of an autoregressive AR(p) model
 
     :math:`x_{k+1}=c+\phi_1 x_k+\phi_2 x_{k-1}+\dots+\phi_p x_{k-p}+\phi_{p+1}\epsilon`
 
@@ -222,7 +222,7 @@ def estimate_ar_params_ols(x, p, d=0, check_stationarity=True,
 def estimate_ar_params_ols_localized(x, p, window_radius, d=0,
                                      include_constant_term=False, h=0, lam=0.0,
                                      window="gaussian"):
-    """Estimate the parameters of a localized AR(p) model
+    r"""Estimate the parameters of a localized AR(p) model
 
     :math:`x_{k+1,i}=c_i+\phi_{1,i}x_{k,i}+\phi_{2,i}x_{k-1,i}+\dots+\phi_{p,i}x_{k-p,i}+\phi_{p+1,i}\epsilon`
 
@@ -365,7 +365,7 @@ def estimate_ar_params_ols_localized(x, p, window_radius, d=0,
 
 
 def estimate_ar_params_yw(gamma, d=0, check_stationarity=True):
-    """Estimate the parameters of an AR(p) model
+    r"""Estimate the parameters of an AR(p) model
 
     :math:`x_{k+1}=\phi_1 x_k+\phi_2 x_{k-1}+\dots+\phi_p x_{k-p}+\phi_{p+1}\epsilon`
 
@@ -441,7 +441,7 @@ def estimate_ar_params_yw(gamma, d=0, check_stationarity=True):
 
 
 def estimate_ar_params_yw_localized(gamma, d=0):
-    """Estimate the parameters of a localized AR(p) model
+    r"""Estimate the parameters of a localized AR(p) model
 
     :math:`x_{k+1,i}=\phi_{1,i}x_{k,i}+\phi_{2,i}x_{k-1,i}+\dots+\phi_{p,i}x_{k-p,i}+\phi_{p+1}\epsilon`
 
@@ -518,7 +518,7 @@ def estimate_ar_params_yw_localized(gamma, d=0):
 
 def estimate_var_params_ols(x, p, d=0, check_stationarity=True,
                             include_constant_term=False, h=0, lam=0.0):
-    """Estimate the parameters of a vector autoregressive VAR(p) model
+    r"""Estimate the parameters of a vector autoregressive VAR(p) model
 
       :math:`\mathbf{x}_{k+1}=\mathbf{c}+\mathbf{\Phi}_1\mathbf{x}_k+
       \mathbf{\Phi}_2\mathbf{x}_{k-1}+\dots+\mathbf{\Phi}_p\mathbf{x}_{k-p}+
@@ -634,7 +634,7 @@ def estimate_var_params_ols(x, p, d=0, check_stationarity=True,
 def estimate_var_params_ols_localized(x, p, window_radius, d=0,
                                       include_constant_term=False,
                                       h=0, lam=0.0, window="gaussian"):
-    """Estimate the parameters of a vector autoregressive VAR(p) model
+    r"""Estimate the parameters of a vector autoregressive VAR(p) model
 
       :math:`\mathbf{x}_{k+1,i}=\mathbf{c}_i+\mathbf{\Phi}_{1,i}\mathbf{x}_{k,i}+
       \mathbf{\Phi}_{2,i}\mathbf{x}_{k-1,i}+\dots+\mathbf{\Phi}_{p,i}
@@ -787,7 +787,7 @@ def estimate_var_params_ols_localized(x, p, window_radius, d=0,
 
 
 def estimate_var_params_yw(gamma, d=0, check_stationarity=True):
-    """Estimate the parameters of a VAR(p) model
+    r"""Estimate the parameters of a VAR(p) model
 
       :math:`\mathbf{x}_{k+1}=\mathbf{\Phi}_1\mathbf{x}_k+
       \mathbf{\Phi}_2\mathbf{x}_{k-1}+\dots+\mathbf{\Phi}_p\mathbf{x}_{k-p}+
@@ -870,7 +870,7 @@ def estimate_var_params_yw(gamma, d=0, check_stationarity=True):
 
 
 def estimate_var_params_yw_localized(gamma, d=0):
-    """Estimate the parameters of a vector autoregressive VAR(p) model
+    r"""Estimate the parameters of a vector autoregressive VAR(p) model
 
       :math:`\mathbf{x}_{k+1,i}=\mathbf{\Phi}_{1,i}\mathbf{x}_{k,i}+
       \mathbf{\Phi}_{2,i}\mathbf{x}_{k-1,i}+\dots+\mathbf{\Phi}_{p,i}
@@ -949,7 +949,7 @@ def estimate_var_params_yw_localized(gamma, d=0):
 
 
 def iterate_ar_model(x, phi, eps=None):
-    """Apply an AR(p) model
+    r"""Apply an AR(p) model
 
     :math:`x_{k+1}=\phi_1 x_k+\phi_2 x_{k-1}+\dots+\phi_p x_{k-p}+\phi_{p+1}\epsilon`
 
@@ -1000,7 +1000,7 @@ def iterate_ar_model(x, phi, eps=None):
 
 
 def iterate_var_model(x, phi, eps=None):
-    """Apply a VAR(p) model
+    r"""Apply a VAR(p) model
 
     :math:`\mathbf{x}_{k+1}=\mathbf{\Phi}_1\mathbf{x}_k+\mathbf{\Phi}_2
     \mathbf{x}_{k-1}+\dots+\mathbf{\Phi}_p\mathbf{x}_{k-p}+
@@ -1057,7 +1057,7 @@ def iterate_var_model(x, phi, eps=None):
 
 
 def test_ar_stationarity(phi):
-    """Test stationarity of an AR(p) process. That is, test that the roots of
+    r"""Test stationarity of an AR(p) process. That is, test that the roots of
     the equation :math:`x^p-\phi_1*x^{p-1}-\dots-\phi_p` lie inside the unit
     circle.
 
@@ -1079,7 +1079,7 @@ def test_ar_stationarity(phi):
 
 
 def test_var_stationarity(phi):
-    """Test stationarity of an AR(p) process. That is, test that the moduli of
+    r"""Test stationarity of an AR(p) process. That is, test that the moduli of
     the eigenvalues of the companion matrix lie inside the unit circle.
 
     Parameters

pysteps/utils/transformation.py

@@ -278,7 +278,7 @@ def NQ_transform(R, metadata=None, inverse=False, **kwargs):
 
     References
     ----------
-    Bogner, K., Pappenberger, field., and Cloke, H. L.: Technical Note: The normal
+    Bogner, K., Pappenberger, F., and Cloke, H. L.: Technical Note: The normal
     quantile transformation and its application in a flood forecasting system,
     Hydrol. Earth Syst. Sci., 16, 1085-1094,
     https://doi.org/10.5194/hess-16-1085-2012, 2012.

pysteps/verification/detcatscores.py

@@ -130,7 +130,7 @@ def det_cat_fct_init(thr, axis=None):
     # catch case of axis passed as integer
     def get_iterable(x):
         if x is None or (
-            isinstance(x, collections.Iterable) and not isinstance(x, int)
+            isinstance(x, collections.abc.Iterable) and not isinstance(x, int)
         ):
             return x
         else:
@@ -332,7 +332,7 @@ def det_cat_fct_compute(contab, scores=""):
 
     # catch case of single score passed as string
     def get_iterable(x):
-        if isinstance(x, collections.Iterable) and not isinstance(x, str):
+        if isinstance(x, collections.abc.Iterable) and not isinstance(x, str):
             return x
         else:
             return (x,)

pysteps/verification/detcontscores.py

@@ -139,7 +139,7 @@ def det_cont_fct(pred, obs, scores="", axis=None, conditioning=None, thr=0.0):
 
     # catch case of single score passed as string
     def get_iterable(x):
-        if isinstance(x, collections.Iterable) and not isinstance(x, str):
+        if isinstance(x, collections.abc.Iterable) and not isinstance(x, str):
             return x
         else:
             return (x,)
@@ -257,7 +257,7 @@ def det_cont_fct_init(axis=None, conditioning=None, thr=0.0):
     # catch case of axis passed as integer
     def get_iterable(x):
         if x is None or (
-            isinstance(x, collections.Iterable) and not isinstance(x, int)
+            isinstance(x, collections.abc.Iterable) and not isinstance(x, int)
         ):
             return x
         else:
@@ -299,7 +299,7 @@ def det_cont_fct_accum(err, pred, obs):
     References
     ----------
 
-    Chan, Tony field.; Golub, Gene H.; LeVeque, Randall J. (1979), "Updating
+    Chan, Tony F.; Golub, Gene H.; LeVeque, Randall J. (1979), "Updating
     Formulae and a Pairwise Algorithm for Computing Sample Variances.",
     Technical Report STAN-CS-79-773, Department of Computer Science,
     Stanford University.
@@ -561,7 +561,7 @@ def det_cont_fct_compute(err, scores=""):
 
     # catch case of single score passed as string
     def get_iterable(x):
-        if isinstance(x, collections.Iterable) and not isinstance(x, str):
+        if isinstance(x, collections.abc.Iterable) and not isinstance(x, str):
             return x
         else:
             return (x,)
@@ -694,7 +694,7 @@ def _scatter(pred, obs, axis=None):
     # catch case of axis passed as integer
     def get_iterable(x):
         if x is None or (
-            isinstance(x, collections.Iterable) and not isinstance(x, int)
+            isinstance(x, collections.abc.Iterable) and not isinstance(x, int)
         ):
             return x
         else:
@@ -748,7 +748,7 @@ def _spearmanr(pred, obs, axis=None):
     # catch case of axis passed as integer
     def get_iterable(x):
         if x is None or (
-            isinstance(x, collections.Iterable) and not isinstance(x, int)
+            isinstance(x, collections.abc.Iterable) and not isinstance(x, int)
         ):
             return x
         else:

pysteps/verification/spatialscores.py

@@ -151,7 +151,7 @@ def intensity_scale_init(name, thrs, scales=None, wavelet="Haar"):
 
     # catch scalars when passed as arguments
     def get_iterable(x):
-        if isinstance(x, collections.Iterable):
+        if isinstance(x, collections.abc.Iterable):
             return np.copy(x)
         else:
             return np.copy((x,))
@@ -453,11 +453,12 @@ def binary_mse_accum(bmse, X_f, X_o):
         mse = np.mean(E_decomp[j] ** 2)
         if np.isfinite(mse):
             bmse["mse"][j] = (bmse["mse"][j] * bmse["n"] + mse) / (
-                bmse["n"] + 1
+                    bmse["n"] + 1
             )
 
     bmse["n"] += 1
 
+
 def binary_mse_merge(bmse_1, bmse_2):
     """Merge two BMSE objects.
 
@@ -503,12 +504,12 @@ def binary_mse_merge(bmse_1, bmse_2):
     # merge the BMSE objects
     bmse = bmse_1.copy()
     bmse["eps"] = (bmse["eps"] * bmse["n"] + bmse_2["eps"] * bmse_2["n"]) / (
-        bmse["n"] + bmse_2["n"]
+            bmse["n"] + bmse_2["n"]
     )
     for j, scale in enumerate(bmse["scales"]):
         bmse["mse"][j] = (
-            bmse["mse"][j] * bmse["n"] + bmse_2["mse"][j] * bmse_2["n"]
-        ) / (bmse["n"] + bmse_2["n"])
+                                 bmse["mse"][j] * bmse["n"] + bmse_2["mse"][j] * bmse_2["n"]
+                         ) / (bmse["n"] + bmse_2["n"])
     bmse["n"] += bmse_2["n"]
 
     return bmse