Cascade scale fix (#283)

* Remove unnecessary encoding lines

* Use default l_0 value that gives constant ratio between scales

* Take l_0 from max(width,height) instead of min

* Fix typo

* Docstring fix

* Use more consistent output variable names

* Add missing output variable description to docstring

* Add option to return callable functions for bandpass filter weights

* Implement alternative strategy for choosing filter weights

* Set the weights corresponding to the first Fourier frequency to zero

* Add option to subtract field mean before decomposition (defaults to True)

* Fix nan values due to division by zero

* Fix typo

* Add the first Fourier wavenumber/field mean to the first filter by default

* Fix typo

* Ensure that weights sum to one

* Set subtract_mean to False by default

* Slightly adjust CRPS thresholds for the new bandpass filter configuration

Co-authored-by: Daniele Nerini <daniele.nerini@gmail.com>

Author: pulkkins

pysteps/cascade/bandpass_filters.py

@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 pysteps.cascade.bandpass_filters
 ================================
@@ -64,10 +63,14 @@ def filter_uniform(shape, n):
     n: int
         Not used. Needed for compatibility with the filter interface.
 
+    Returns
+    -------
+    out: dict
+        A dictionary containing the filter.
     """
     del n  # Unused
 
-    result = {}
+    out = {}
 
     try:
         height, width = shape
@@ -76,17 +79,23 @@ def filter_uniform(shape, n):
 
     r_max = int(max(width, height) / 2) + 1
 
-    result["weights_1d"] = np.ones((1, r_max))
-    result["weights_2d"] = np.ones((1, height, int(width / 2) + 1))
-    result["central_freqs"] = None
-    result["central_wavenumbers"] = None
-    result["shape"] = shape
+    out["weights_1d"] = np.ones((1, r_max))
+    out["weights_2d"] = np.ones((1, height, int(width / 2) + 1))
+    out["central_freqs"] = None
+    out["central_wavenumbers"] = None
+    out["shape"] = shape
 
-    return result
+    return out
 
 
 def filter_gaussian(
-    shape, n, l_0=3, gauss_scale=0.5, gauss_scale_0=0.5, d=1.0, normalize=True
+    shape,
+    n,
+    gauss_scale=0.5,
+    d=1.0,
+    normalize=True,
+    return_weight_funcs=False,
+    include_mean=True,
 ):
     """
     Implements a set of Gaussian bandpass filters in logarithmic frequency
@@ -99,20 +108,20 @@ def filter_gaussian(
         the domain is assumed to have square shape.
     n: int
         The number of frequency bands to use. Must be greater than 2.
-    l_0: int
-        Central frequency of the second band (the first band is always centered
-        at zero).
     gauss_scale: float
-        Optional scaling prameter. Proportional to the standard deviation of
+        Optional scaling parameter. Proportional to the standard deviation of
         the Gaussian weight functions.
-    gauss_scale_0: float
-        Optional scaling parameter for the Gaussian function corresponding to
-        the first frequency band.
     d: scalar, optional
         Sample spacing (inverse of the sampling rate). Defaults to 1.
     normalize: bool
         If True, normalize the weights so that for any given wavenumber
         they sum to one.
+    return_weight_funcs: bool
+        If True, add callable weight functions to the output dictionary with
+        the key 'weight_funcs'.
+    include_mean: bool
+        If True, include the first Fourier wavenumber (corresponding to the
+        field mean) to the first filter.
 
     Returns
     -------
@@ -133,6 +142,8 @@ def filter_gaussian(
     except TypeError:
         height, width = (shape, shape)
 
+    max_length = max(width, height)
+
     rx = np.s_[: int(width / 2) + 1]
 
     if (height % 2) == 1:
@@ -145,13 +156,13 @@ def filter_gaussian(
 
     r_2d = np.roll(np.sqrt(x_grid * x_grid + y_grid * y_grid), dy, axis=0)
 
-    max_length = max(width, height)
-
     r_max = int(max_length / 2) + 1
     r_1d = np.arange(r_max)
 
     wfs, central_wavenumbers = _gaussweights_1d(
-        max_length, n, l_0=l_0, gauss_scale=gauss_scale, gauss_scale_0=gauss_scale_0
+        max_length,
+        n,
+        gauss_scale=gauss_scale,
     )
 
     weights_1d = np.empty((n, r_max))
@@ -168,36 +179,48 @@ def filter_gaussian(
             weights_1d[k, :] /= weights_1d_sum
             weights_2d[k, :, :] /= weights_2d_sum
 
-    result = {"weights_1d": weights_1d, "weights_2d": weights_2d}
-    result["shape"] = shape
+    for i in range(len(wfs)):
+        if i == 0 and include_mean:
+            weights_1d[i, 0] = 1.0
+            weights_2d[i, 0, 0] = 1.0
+        else:
+            weights_1d[i, 0] = 0.0
+            weights_2d[i, 0, 0] = 0.0
+
+    out = {"weights_1d": weights_1d, "weights_2d": weights_2d}
+    out["shape"] = shape
 
     central_wavenumbers = np.array(central_wavenumbers)
-    result["central_wavenumbers"] = central_wavenumbers
+    out["central_wavenumbers"] = central_wavenumbers
 
     # Compute frequencies
     central_freqs = 1.0 * central_wavenumbers / max_length
     central_freqs[0] = 1.0 / max_length
     central_freqs[-1] = 0.5  # Nyquist freq
     central_freqs = 1.0 * d * central_freqs
-    result["central_freqs"] = central_freqs
+    out["central_freqs"] = central_freqs
+
+    if return_weight_funcs:
+        out["weight_funcs"] = wfs
 
-    return result
+    return out
 
 
-def _gaussweights_1d(l, n, l_0=3, gauss_scale=0.5, gauss_scale_0=0.5):
-    e = pow(0.5 * l / l_0, 1.0 / (n - 2))
-    r = [(l_0 * pow(e, k - 1), l_0 * pow(e, k)) for k in range(1, n - 1)]
+def _gaussweights_1d(l, n, gauss_scale=0.5):
+    q = pow(0.5 * l, 1.0 / n)
+    r = [(pow(q, k - 1), pow(q, k)) for k in range(1, n + 1)]
+    r = [0.5 * (r_[0] + r_[1]) for r_ in r]
 
     def log_e(x):
         if len(np.shape(x)) > 0:
             res = np.empty(x.shape)
             res[x == 0] = 0.0
-            res[x > 0] = np.log(x[x > 0]) / np.log(e)
+            res[x > 0] = np.log(x[x > 0]) / np.log(q)
         else:
             if x == 0.0:
                 res = 0.0
             else:
-                res = np.log(x) / np.log(e)
+                res = np.log(x) / np.log(q)
 
         return res
 
@@ -211,25 +234,11 @@ def __call__(self, x):
             return np.exp(-(x**2.0) / (2.0 * self.s**2.0))
 
     weight_funcs = []
-    central_wavenumbers = [0.0]
-
-    weight_funcs.append(GaussFunc(0.0, gauss_scale_0))
+    central_wavenumbers = []
 
     for i, ri in enumerate(r):
-        rc = log_e(ri[0])
+        rc = log_e(ri)
         weight_funcs.append(GaussFunc(rc, gauss_scale))
-        central_wavenumbers.append(ri[0])
-
-    gf = GaussFunc(log_e(l / 2), gauss_scale)
-
-    def g(x):
-        res = np.ones(x.shape)
-        mask = x <= l / 2
-        res[mask] = gf(x[mask])
-
-        return res
-
-    weight_funcs.append(g)
-    central_wavenumbers.append(l / 2)
+        central_wavenumbers.append(ri)
 
     return weight_funcs, central_wavenumbers

pysteps/cascade/decomposition.py

@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 pysteps.cascade.decomposition
 =============================
@@ -14,9 +13,8 @@
 where field is the input field and bp_filter is a dictionary returned by a
 filter method implemented in :py:mod:`pysteps.cascade.bandpass_filters`. The
 decomp argument is a decomposition obtained by calling decomposition_xxx.
-Optional parameters can be passed in
-the keyword arguments. The output of each method is a dictionary with the
-following key-value pairs:
+Optional parameters can be passed in the keyword arguments. The output of each
+method is a dictionary with the following key-value pairs:
 
 +-------------------+----------------------------------------------------------+
 |        Key        |                      Value                               |
@@ -120,6 +118,10 @@ def decomposition_fft(field, bp_filter, **kwargs):
         Applicable if output_domain is "spectral". If set to True, only the
         parts of the Fourier spectrum with non-negligible filter weights are
         stored. Defaults to False.
+    subtract_mean: bool
+        If set to True, subtract the mean value before the decomposition and
+        store it to the output dictionary. Applicable if input_domain is
+        "spatial". Defaults to False.
 
     Returns
     -------
@@ -138,6 +140,7 @@ def decomposition_fft(field, bp_filter, **kwargs):
     output_domain = kwargs.get("output_domain", "spatial")
     compute_stats = kwargs.get("compute_stats", True)
     compact_output = kwargs.get("compact_output", False)
+    subtract_mean = kwargs.get("subtract_mean", False)
 
     if normalize and not compute_stats:
         compute_stats = True
@@ -194,6 +197,11 @@ def decomposition_fft(field, bp_filter, **kwargs):
     means = []
     stds = []
 
+    if subtract_mean and input_domain == "spatial":
+        field_mean = np.mean(field)
+        field = field - field_mean
+        result["field_mean"] = field_mean
+
     if input_domain == "spatial":
         field_fft = fft.rfft2(field)
     else:
@@ -276,7 +284,7 @@ def recompose_fft(decomp, **kwargs):
     if not decomp["normalized"] and not (
         decomp["domain"] == "spectral" and decomp["compact_output"]
     ):
-        return np.sum(levels, axis=0)
+        result = np.sum(levels, axis=0)
     else:
         if decomp["compact_output"]:
             weight_masks = decomp["weight_masks"]
@@ -291,4 +299,7 @@ def recompose_fft(decomp, **kwargs):
             result = [levels[i] * sigma[i] + mu[i] for i in range(len(levels))]
             result = np.sum(np.stack(result), axis=0)
 
-        return result
+    if "field_mean" in decomp:
+        result += decomp["field_mean"]
+
+    return result

pysteps/cascade/interface.py

@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 pysteps.cascade.interface
 =========================

pysteps/tests/test_nowcasts_steps.py

@@ -24,9 +24,9 @@
 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.25),
-    (5, 6, 2, "sprog", None, "spatial", 3, 8.35),
-    (5, 6, 2, "obs", None, "spatial", 3, 8.30),
+    (5, 6, 2, "incremental", None, "spatial", 3, 7.31),
+    (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),
     (5, 6, 2, None, "mean", "spatial", 3, 1.35),
     (5, 6, 2, "incremental", "cdf", "spectral", 3, 0.60),