Improve documentation of rainfarm downscaling (#292)

Author: dnerini

examples/rainfarm_downscale.py

@@ -57,11 +57,15 @@
 # -----------------
 #
 # To test our downscaling method, we first need to upscale the original field to
-# a lower resolution. We are going to use an upscaling factor of 16 x.
+# a lower resolution. This is only for demo purposes, as we need to artificially
+# create a lower resolution field to apply our downscaling method.
+# We are going to use a factor of 16 x.
 
-upscaling_factor = 16
-upscale_to = metadata["xpixelsize"] * upscaling_factor  # upscaling factor : 16 x
-precip_lr, metadata_lr = aggregate_fields_space(precip, metadata, upscale_to)
+scale_factor = 16
+upscaled_resolution = (
+    metadata["xpixelsize"] * scale_factor
+)  # upscaled resolution : 16 km
+precip_lr, metadata_lr = aggregate_fields_space(precip, metadata, upscaled_resolution)
 
 # Plot the upscaled rainfall field
 plt.figure()
@@ -81,33 +85,29 @@
 # Per realization, generate a stochastically downscaled precipitation field
 # and plot it.
 # The first time, the spectral slope alpha needs to be estimated. To illustrate
-# the sensitity of this parameter, we are going to plot some realizations with
+# the sensitivity of this parameter, we are going to plot some realizations with
 # half or double the estimated slope.
 alpha = None
 for n in range(num_realizations):
 
     # Spectral slope estimated from the upscaled field
     precip_hr, alpha = rainfarm.downscale(
-        precip_lr, alpha=alpha, ds_factor=upscaling_factor, return_alpha=True
+        precip_lr, ds_factor=scale_factor, alpha=alpha, return_alpha=True
     )
     plt.subplot(num_realizations, 3, n * 3 + 2)
     plot_precip_field(precip_hr, geodata=metadata, axis="off", colorbar=False)
     if n == 0:
         plt.title(f"alpha={alpha:.1f}")
 
     # Half the estimated slope
-    precip_hr = rainfarm.downscale(
-        precip_lr, alpha=alpha * 0.5, ds_factor=upscaling_factor
-    )
+    precip_hr = rainfarm.downscale(precip_lr, ds_factor=scale_factor, alpha=alpha * 0.5)
     plt.subplot(num_realizations, 3, n * 3 + 1)
     plot_precip_field(precip_hr, geodata=metadata, axis="off", colorbar=False)
     if n == 0:
         plt.title(f"alpha={alpha * 0.5:.1f}")
 
     # Double the estimated slope
-    precip_hr = rainfarm.downscale(
-        precip_lr, alpha=alpha * 2, ds_factor=upscaling_factor
-    )
+    precip_hr = rainfarm.downscale(precip_lr, ds_factor=scale_factor, alpha=alpha * 2)
     plt.subplot(num_realizations, 3, n * 3 + 3)
     plot_precip_field(precip_hr, geodata=metadata, axis="off", colorbar=False)
     if n == 0:

pysteps/downscaling/rainfarm.py

@@ -6,6 +6,11 @@
 Implementation of the RainFARM stochastic downscaling method as described in
 :cite:`Rebora2006`.
 
+RainFARM is a downscaling algorithm for rainfall fields developed by Rebora et
+al. (2006). The method can represent the realistic small-scale variability of the
+downscaled precipitation field by means of Gaussian random fields.
+
+
 .. autosummary::
     :toctree: ../generated/
 
@@ -39,38 +44,35 @@ def _balanced_spatial_average(x, k):
     return convolve(x, k) / convolve(ones, k)
 
 
-def downscale(precip, alpha=None, ds_factor=16, threshold=None, return_alpha=False):
+def downscale(precip, ds_factor, alpha=None, threshold=None, return_alpha=False):
     """
-    Downscale a rainfall field by a given factor.
+    Downscale a rainfall field by increasing its spatial resolution by
+    a positive integer factor.
 
     Parameters
     ----------
-
-    precip: array_like
-        Array of shape (m,n) containing the input field.
+    precip: array-like
+        Array of shape (m, n) containing the input field.
         The input is expected to contain rain rate values.
-
+        All values are required to be finite.
     alpha: float, optional
-        Spectral slope. If none, the slope is estimated from
+        Spectral slope. If None, the slope is estimated from
         the input array.
-
-    ds_factor: int, optional
-        Downscaling factor.
-
+    ds_factor: positive int
+        Downscaling factor, it specifies by how many times
+        to increase the initial grid resolution.
     threshold: float, optional
         Set all values lower than the threshold to zero.
-
     return_alpha: bool, optional
-        Whether to return the estimated spectral slope `alpha`.
-
+        Whether to return the estimated spectral slope ``alpha``.
 
     Returns
     -------
-    r: array_like
-        Array of shape (m*ds_factor,n*ds_factor) containing
+    r: array-like
+        Array of shape (m * ds_factor, n * ds_factor) containing
         the downscaled field.
     alpha: float
-        Returned only when `return_alpha=True`.
+        Returned only when ``return_alpha=True``.
 
     Notes
     -----

pysteps/tests/test_downscaling_rainfarm.py

@@ -7,30 +7,35 @@
 from pysteps.utils import aggregate_fields_space, square_domain
 
 
-# load and preprocess input field
-precip, metadata = get_precipitation_fields(
-    num_prev_files=0, num_next_files=0, return_raw=False, metadata=True
-)
-precip = precip.filled()
-precip, metadata = square_domain(precip, metadata, "crop")
+@pytest.fixture(scope="module")
+def data():
+    precip, metadata = get_precipitation_fields(
+        num_prev_files=0, num_next_files=0, return_raw=False, metadata=True
+    )
+    precip = precip.filled()
+    precip, metadata = square_domain(precip, metadata, "crop")
+    return precip, metadata
 
 
 rainfarm_arg_names = ("alpha", "ds_factor", "threshold", "return_alpha")
-
-
 rainfarm_arg_values = [(1.0, 1, 0, False), (1, 2, 0, False), (1, 4, 0, False)]
 
 
 @pytest.mark.parametrize(rainfarm_arg_names, rainfarm_arg_values)
-def test_rainfarm_shape(alpha, ds_factor, threshold, return_alpha):
-    """Test that the output of rainfarm is consistent with the downscalnig factor."""
-
+def test_rainfarm_shape(data, alpha, ds_factor, threshold, return_alpha):
+    """Test that the output of rainfarm is consistent with the downscaling factor."""
+    precip, metadata = data
     window = metadata["xpixelsize"] * ds_factor
     precip_lr, __ = aggregate_fields_space(precip, metadata, window)
 
     rainfarm = downscaling.get_method("rainfarm")
-
-    precip_hr = rainfarm(precip_lr, alpha, ds_factor, threshold, return_alpha)
+    precip_hr = rainfarm(
+        precip_lr,
+        alpha=alpha,
+        ds_factor=ds_factor,
+        threshold=threshold,
+        return_alpha=return_alpha,
+    )
 
     assert precip_hr.ndim == precip.ndim
     assert precip_hr.shape[0] == precip.shape[0]
@@ -41,15 +46,20 @@ def test_rainfarm_shape(alpha, ds_factor, threshold, return_alpha):
 
 
 @pytest.mark.parametrize(rainfarm_arg_names, rainfarm_arg_values)
-def test_rainfarm_alpha(alpha, ds_factor, threshold, return_alpha):
+def test_rainfarm_alpha(data, alpha, ds_factor, threshold, return_alpha):
     """Test that rainfarm computes and returns alpha."""
-
+    precip, metadata = data
     window = metadata["xpixelsize"] * ds_factor
     precip_lr, __ = aggregate_fields_space(precip, metadata, window)
 
     rainfarm = downscaling.get_method("rainfarm")
-
-    precip_hr = rainfarm(precip_lr, alpha, ds_factor, threshold, return_alpha)
+    precip_hr = rainfarm(
+        precip_lr,
+        alpha=alpha,
+        ds_factor=ds_factor,
+        threshold=threshold,
+        return_alpha=return_alpha,
+    )
 
     assert len(precip_hr) == 2
     if alpha is None: