Skip to content

Commit

Permalink
Improve documentation of rainfarm downscaling (#292)
Browse files Browse the repository at this point in the history
  • Loading branch information
dnerini authored Aug 6, 2022
1 parent 33d80d4 commit ed76516
Show file tree
Hide file tree
Showing 3 changed files with 58 additions and 46 deletions.
24 changes: 12 additions & 12 deletions examples/rainfarm_downscale.py
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand All @@ -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:
Expand Down
36 changes: 19 additions & 17 deletions pysteps/downscaling/rainfarm.py
Original file line number Diff line number Diff line change
Expand Up @@ -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/
Expand Down Expand Up @@ -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
-----
Expand Down
44 changes: 27 additions & 17 deletions pysteps/tests/test_downscaling_rainfarm.py
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand All @@ -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:
Expand Down

0 comments on commit ed76516

Please sign in to comment.