Update RainFARM example following #311 (#367)

* Update RainFARM example following #311

examples/rainfarm_downscale.py

@@ -10,43 +10,71 @@
 al. (2006). The method can represent the realistic small-scale variability of the
 downscaled precipitation field by means of Gaussian random fields.
 
+Steps:
+ 1. Read the input precipitation data.
+ 2. Upscale the precipitation field.
+ 3. Downscale the field to its original resolution using RainFARM with defaults.
+ 4. Downscale with smoothing.
+ 5. Downscale with spectral fusion.
+ 6. Downscale with smoothing and spectral fusion.
+
+References:
+
+ Rebora, N., L. Ferraris, J. von Hardenberg, and A. Provenzale, 2006: RainFARM:
+ Rainfall downscaling by a filtered autoregressive model. J. Hydrometeor., 7,
+ 724–738.
+
+ D D'Onofrio, E Palazzi, J von Hardenberg, A Provenzale, and S Calmanti, 2014:
+ Stochastic rainfall downscaling of climate models. J. Hydrometeorol., 15(2):830–843.
 """
 
 import matplotlib.pyplot as plt
 import numpy as np
 import os
 from pprint import pprint
+import logging
 
 from pysteps import io, rcparams
 from pysteps.utils import aggregate_fields_space, square_domain, to_rainrate
 from pysteps.downscaling import rainfarm
 from pysteps.visualization import plot_precip_field
 
+# Configure logging
+logging.basicConfig(
+ level=logging.INFO, format="%(asctime)s - %(levelname)s - %(message)s"
+)
+
 ###############################################################################
 # Read the input data
 # -------------------
 #
 # As first step, we need to import the precipitation field that we are going
 # to use in this example.
 
+
+def read_precipitation_data(file_path):
+ """Read and process precipitation data from a file."""
+ precip, _, metadata = io.import_mch_gif(
+ file_path, product="AQC", unit="mm", accutime=5.0
+ )
+ precip, metadata = to_rainrate(precip, metadata)
+ precip, metadata = square_domain(precip, metadata, "crop")
+ return precip, metadata
+
+
 # Import the example radar composite
 root_path = rcparams.data_sources["mch"]["root_path"]
 filename = os.path.join(root_path, "20160711", "AQC161932100V_00005.801.gif")
-precip, _, metadata = io.import_mch_gif(
- filename, product="AQC", unit="mm", accutime=5.0
-)
-
-# Convert to mm/h
-precip, metadata = to_rainrate(precip, metadata)
 
-# Reduce to a square domain
-precip, metadata = square_domain(precip, metadata, "crop")
+# Read and process data
+precip, metadata = read_precipitation_data(filename)
 
 # Nicely print the metadata
 pprint(metadata)
 
 # Plot the original rainfall field
 plot_precip_field(precip, geodata=metadata)
+plt.title("Original Rainfall Field")
 plt.show()
 
 # Assign the fill value to all the Nans
@@ -61,60 +89,87 @@
 # create a lower resolution field to apply our downscaling method.
 # We are going to use a factor of 16 x.
 
+
+def upscale_field(precip, metadata, scale_factor):
+ """Upscale the precipitation field by a given scale factor."""
+ upscaled_resolution = metadata["xpixelsize"] * scale_factor
+ precip_lr, metadata_lr = aggregate_fields_space(
+ precip, metadata, upscaled_resolution
+ )
+ return precip_lr, metadata_lr
+
+
 scale_factor = 16
-upscaled_resolution = (
- metadata["xpixelsize"] * scale_factor
-) # upscaled resolution : 16 km
-precip_lr, metadata_lr = aggregate_fields_space(precip, metadata, upscaled_resolution)
+precip_lr, metadata_lr = upscale_field(precip, metadata, scale_factor)
 
 # Plot the upscaled rainfall field
 plt.figure()
 plot_precip_field(precip_lr, geodata=metadata_lr)
+plt.title("Upscaled Rainfall Field")
+plt.show()
 
 ###############################################################################
 # Downscale the field
 # -------------------
 #
-# We can now use RainFARM to generate stochastic realizations of the downscaled
-# precipitation field.
-
-fig = plt.figure(figsize=(5, 8))
-# Set the number of stochastic realizations
-num_realizations = 5
-
-# 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 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, 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, 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, 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:
- plt.title(f"alpha={alpha * 2:.1f}")
-
- plt.subplots_adjust(wspace=0, hspace=0)
-
-plt.tight_layout()
+# We can now use RainFARM to downscale the precipitation field.
+
+# Basic downscaling
+precip_hr = rainfarm.downscale(precip_lr, ds_factor=scale_factor)
+
+# Plot the downscaled rainfall field
+plt.figure()
+plot_precip_field(precip_hr, geodata=metadata)
+plt.title("Downscaled Rainfall Field")
+plt.show()
+
+###############################################################################
+# Downscale with smoothing
+# ------------------------
+#
+# Add smoothing with a Gaussian kernel during the downscaling process.
+
+precip_hr_smooth = rainfarm.downscale(
+ precip_lr, ds_factor=scale_factor, kernel_type="gaussian"
+)
+
+# Plot the downscaled rainfall field with smoothing
+plt.figure()
+plot_precip_field(precip_hr_smooth, geodata=metadata)
+plt.title("Downscaled Rainfall Field with Gaussian Smoothing")
+plt.show()
+
+###############################################################################
+# Downscale with spectral fusion
+# ------------------------------
+#
+# Apply spectral merging as described in D'Onofrio et al. (2014).
+
+precip_hr_fusion = rainfarm.downscale(
+ precip_lr, ds_factor=scale_factor, spectral_fusion=True
+)
+
+# Plot the downscaled rainfall field with spectral fusion
+plt.figure()
+plot_precip_field(precip_hr_fusion, geodata=metadata)
+plt.title("Downscaled Rainfall Field with Spectral Fusion")
+plt.show()
+
+###############################################################################
+# Combined Downscale with smoothing and spectral fusion
+# -----------------------------------------------------
+#
+# Apply both smoothing with a Gaussian kernel and spectral fusion during the
+# downscaling process to observe the combined effect.
+
+precip_hr_combined = rainfarm.downscale(
+ precip_lr, ds_factor=scale_factor, kernel_type="gaussian", spectral_fusion=True
+)
+
+# Plot the downscaled rainfall field with smoothing and spectral fusion
+plt.figure()
+plot_precip_field(precip_hr_combined, geodata=metadata)
+plt.title("Downscaled Rainfall Field with Gaussian Smoothing and Spectral Fusion")
 plt.show()
 
 ###############################################################################
@@ -126,13 +181,4 @@
 # the original algorithm from Rebora et al. (2006), it cannot downscale the temporal
 # dimension.
 
-
-###############################################################################
-# References
-# ----------
-#
-# Rebora, N., L. Ferraris, J. von Hardenberg, and A. Provenzale, 2006: RainFARM:
-# Rainfall downscaling by a filtered autoregressive model. J. Hydrometeor., 7,
-# 724–738.
-
 # sphinx_gallery_thumbnail_number = 2