From 787e39c1d82634d2f76c616f9b9d28722c4a84fb Mon Sep 17 00:00:00 2001 From: Simon Beylat <72618409+simonbeylat@users.noreply.github.com> Date: Wed, 21 Sep 2022 14:06:04 +0200 Subject: [PATCH 01/21] Update rainfarm.py Add smoothing operator --- pysteps/downscaling/rainfarm.py | 70 +++++++++++++++++++++++++++++++-- 1 file changed, 67 insertions(+), 3 deletions(-) diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index 4b9bb0bf7..151795dad 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -44,7 +44,60 @@ def _balanced_spatial_average(x, k): return convolve(x, k) / convolve(ones, k) -def downscale(precip, ds_factor, alpha=None, threshold=None, return_alpha=False): +def _smoothconv(P,nas): + + """ + Parameters + ---------- + + P: matrix + matrix with the input field to smoothen, with dimensions ns*ns + + nas : int + original size + + Returns + ------- + + The smoothened field. + + References + ---------- + + Terzago et al. 2018 + + """ + + print("smoothing") + indmask= ~ np.isfinite(P) + P[indmask] = 0. + ns = np.shape(P)[1] + + sdim = (ns/nas)/2 + mask=np.zeros([ns,ns]) + for i in range(ns): + for j in range(ns): + + kx = i + ky = j + if i > (ns/2): + kx=i-ns + if j > (ns/2): + ky=j-ns + r2 = kx * kx + ky * ky + mask[i,j] = np.exp(-(r2 / (sdim*sdim))/2) + + fm = np.fft.fft2(mask) + pf = np.real(np.fft.ifft2(fm*np.fft.fft2(P))) / np.sum(mask) + if np.sum(indmask) > 0: + P[~ indmask] = 1 + pf = pf / (np.real(np.fft.ifft2(fm*np.fft.fft2(P))) / np.sum(P) / len(fm)) + + pf[indmask]=np.nan + return pf + + +def downscale(precip, ds_factor, alpha=None, threshold=None, return_alpha=False , smooth=False): """ Downscale a rainfall field by increasing its spatial resolution by a positive integer factor. @@ -86,6 +139,7 @@ def downscale(precip, ds_factor, alpha=None, threshold=None, return_alpha=False) :cite:`Rebora2006` """ + nas = np.shape(precip)[1] ki = np.fft.fftfreq(precip.shape[0]) kj = np.fft.fftfreq(precip.shape[1]) @@ -119,8 +173,18 @@ def downscale(precip, ds_factor, alpha=None, threshold=None, return_alpha=False) tophat = ((mx**2 + my**2) <= rad**2).astype(float) tophat /= tophat.sum() - P_agg = _balanced_spatial_average(P_u, tophat) - r_agg = _balanced_spatial_average(r, tophat) + if smooth: + P_agg = _balanced_spatial_average(P_u, tophat) + r_agg = _balanced_spatial_average(r, tophat) + + P_agg = _smoothconv(P_agg, nas) + r_agg = _smoothconv(r_agg, nas) + + else: + P_agg = _balanced_spatial_average(P_u, tophat) + r_agg = _balanced_spatial_average(r, tophat) + + r *= P_agg / r_agg if threshold is not None: From e10c7ce2f31620a0426d09401e19948079263446 Mon Sep 17 00:00:00 2001 From: simonbeylat Date: Mon, 3 Oct 2022 16:17:40 +0200 Subject: [PATCH 02/21] add smoothing option {None,'Gaussian','tophat'}, update the _balanced_spatial_average method to add the treatment of missing values, add None option, black autoformatter, add Terzago ref --- doc/source/references.bib | 12 ++++ pysteps/downscaling/rainfarm.py | 117 ++++++++++++++++++++------------ 2 files changed, 85 insertions(+), 44 deletions(-) diff --git a/doc/source/references.bib b/doc/source/references.bib index 66b4b9720..9f5eaf532 100644 --- a/doc/source/references.bib +++ b/doc/source/references.bib @@ -355,3 +355,15 @@ @ARTICLE{ZR2009 YEAR = 2009, DOI = "10.1016/j.atmosres.2009.03.004" } + + +@Article{Terzago2018, +AUTHOR = "Terzago, S. and Palazzi, E. and von Hardenberg, J.", +TITLE = "Stochastic downscaling of precipitation in complex orography: a simple method to reproduce a realistic fine-scale climatology", +JOURNAL = "Natural Hazards and Earth System Sciences", +VOLUME = 18, +YEAR = 2018, +NUMBER = 11, +PAGES = "2825--2840", +DOI = "10.5194/nhess-18-2825-2018" +} \ No newline at end of file diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index 151795dad..36ea9ba0e 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -21,6 +21,7 @@ import numpy as np from scipy.ndimage import convolve +from scipy.ndimage import zoom def _log_slope(log_k, log_power_spectrum): @@ -40,64 +41,82 @@ def _log_slope(log_k, log_power_spectrum): def _balanced_spatial_average(x, k): + indmask = ~np.isfinite(x) + x[indmask] = 0.0 + ones = np.ones_like(x) return convolve(x, k) / convolve(ones, k) -def _smoothconv(P,nas): +def _smoothconv(pr_high_res, orig_res): """ Parameters ---------- - - P: matrix - matrix with the input field to smoothen, with dimensions ns*ns - - nas : int - original size + pr_high_res: matrix + matrix with the input field to smoothen, with dimensions ns*ns. + orig_res : int + original size of the precipitation field. Returns ------- - The smoothened field. References ---------- - - Terzago et al. 2018 + :cite:`Terzago2018` """ - print("smoothing") - indmask= ~ np.isfinite(P) - P[indmask] = 0. - ns = np.shape(P)[1] + indmask = ~np.isfinite(pr_high_res) + pr_high_res[indmask] = 0.0 + ns = np.shape(pr_high_res)[1] - sdim = (ns/nas)/2 - mask=np.zeros([ns,ns]) - for i in range(ns): + sdim = (ns / orig_res) / 2 + mask = np.zeros([ns, ns]) + for i in range(ns): for j in range(ns): - + kx = i ky = j - if i > (ns/2): - kx=i-ns - if j > (ns/2): - ky=j-ns - r2 = kx * kx + ky * ky - mask[i,j] = np.exp(-(r2 / (sdim*sdim))/2) + if i > (ns / 2): + kx = i - ns + if j > (ns / 2): + ky = j - ns + r2 = kx * kx + ky * ky + mask[i, j] = np.exp(-(r2 / (sdim * sdim)) / 2) fm = np.fft.fft2(mask) - pf = np.real(np.fft.ifft2(fm*np.fft.fft2(P))) / np.sum(mask) - if np.sum(indmask) > 0: - P[~ indmask] = 1 - pf = pf / (np.real(np.fft.ifft2(fm*np.fft.fft2(P))) / np.sum(P) / len(fm)) - - pf[indmask]=np.nan - return pf + pf = np.real(np.fft.ifft2(fm * np.fft.fft2(pr_high_res))) / np.sum(mask) + if np.sum(indmask) > 0: + pr_high_res[~indmask] = 1 + pf = pf / ( + np.real(np.fft.ifft2(fm * np.fft.fft2(pr_high_res))) + / np.sum(pr_high_res) + / len(fm) + ) + + pf[indmask] = np.nan + return pf + + +def _check_smooth_value(smooth): + + if smooth is None: + return 0 + elif smooth == "Gaussian": + return 1 + elif smooth == "tophat": + return 2 + else: + raise RuntimeError( + 'the smooth value does not exist, choose from this list {None,"Gaussian","tophat"}' + ) -def downscale(precip, ds_factor, alpha=None, threshold=None, return_alpha=False , smooth=False): +def downscale( + precip, ds_factor, alpha=None, threshold=None, return_alpha=False, smooth=None +): """ Downscale a rainfall field by increasing its spatial resolution by a positive integer factor. @@ -118,7 +137,8 @@ def downscale(precip, ds_factor, alpha=None, threshold=None, return_alpha=False Set all values lower than the threshold to zero. return_alpha: bool, optional Whether to return the estimated spectral slope ``alpha``. - + smooth: list, optional + Add the smoothing operatore from this list {None,"Gaussian","tophat"} Returns ------- r: array-like @@ -139,7 +159,9 @@ def downscale(precip, ds_factor, alpha=None, threshold=None, return_alpha=False :cite:`Rebora2006` """ - nas = np.shape(precip)[1] + type_smooth = _check_smooth_value(smooth) + + orig_res = np.shape(precip)[1] ki = np.fft.fftfreq(precip.shape[0]) kj = np.fft.fftfreq(precip.shape[1]) @@ -173,19 +195,26 @@ def downscale(precip, ds_factor, alpha=None, threshold=None, return_alpha=False tophat = ((mx**2 + my**2) <= rad**2).astype(float) tophat /= tophat.sum() - if smooth: - P_agg = _balanced_spatial_average(P_u, tophat) - r_agg = _balanced_spatial_average(r, tophat) - - P_agg = _smoothconv(P_agg, nas) - r_agg = _smoothconv(r_agg, nas) + if type_smooth == 0: - else: - P_agg = _balanced_spatial_average(P_u, tophat) - r_agg = _balanced_spatial_average(r, tophat) + P_agg = precip + r_agg = zoom(r, 1 / ds_factor, order=1) + factor = np.repeat( + np.repeat(P_agg / r_agg, ds_factor, axis=0), ds_factor, axis=1 + ) + r *= factor + elif type_smooth == 1: - r *= P_agg / r_agg + P_agg = _smoothconv(P_u, orig_res) + r_agg = _smoothconv(r, orig_res) + r *= P_agg / r_agg + + elif type_smooth == 2: + + P_agg = _balanced_spatial_average(P_u, tophat) + r_agg = _balanced_spatial_average(r, tophat) + r *= P_agg / r_agg if threshold is not None: r[r < threshold] = 0 From b67f8c9d3be472c6e7c4d6c1df5d44ef1b33a95e Mon Sep 17 00:00:00 2001 From: simonbeylat Date: Fri, 14 Oct 2022 17:17:07 +0200 Subject: [PATCH 03/21] merge _smoothconv and _balanced_spatial_average function --- pysteps/downscaling/rainfarm.py | 88 ++++++++++++++++++--------------- 1 file changed, 47 insertions(+), 41 deletions(-) diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index 36ea9ba0e..e6f83b0f3 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -40,23 +40,19 @@ def _log_slope(log_k, log_power_spectrum): return alpha -def _balanced_spatial_average(x, k): - indmask = ~np.isfinite(x) - x[indmask] = 0.0 - - ones = np.ones_like(x) - return convolve(x, k) / convolve(ones, k) - - -def _smoothconv(pr_high_res, orig_res): +def _balanced_spatial_average(pr_high_res, tophat=None, orig_res=None, kernel="tophat"): """ Parameters ---------- - pr_high_res: matrix + pr_high_res : matrix matrix with the input field to smoothen, with dimensions ns*ns. + tophat : array_like + Array of weights orig_res : int original size of the precipitation field. + kernel : string + Choose the smoothing method Returns ------- @@ -70,34 +66,44 @@ def _smoothconv(pr_high_res, orig_res): indmask = ~np.isfinite(pr_high_res) pr_high_res[indmask] = 0.0 - ns = np.shape(pr_high_res)[1] - - sdim = (ns / orig_res) / 2 - mask = np.zeros([ns, ns]) - for i in range(ns): - for j in range(ns): - - kx = i - ky = j - if i > (ns / 2): - kx = i - ns - if j > (ns / 2): - ky = j - ns - r2 = kx * kx + ky * ky - mask[i, j] = np.exp(-(r2 / (sdim * sdim)) / 2) - - fm = np.fft.fft2(mask) - pf = np.real(np.fft.ifft2(fm * np.fft.fft2(pr_high_res))) / np.sum(mask) - if np.sum(indmask) > 0: - pr_high_res[~indmask] = 1 - pf = pf / ( - np.real(np.fft.ifft2(fm * np.fft.fft2(pr_high_res))) - / np.sum(pr_high_res) - / len(fm) - ) - pf[indmask] = np.nan - return pf + if kernel == "tophat": + ones = np.ones_like(pr_high_res) + return convolve(pr_high_res, tophat) / convolve(ones, tophat) + + elif kernel == "Gaussian": + ns = np.shape(pr_high_res)[1] + sdim = (ns / orig_res) / 2 + mask = np.zeros([ns, ns]) + for i in range(ns): + for j in range(ns): + + kx = i + ky = j + if i > (ns / 2): + kx = i - ns + if j > (ns / 2): + ky = j - ns + r2 = kx * kx + ky * ky + mask[i, j] = np.exp(-(r2 / (sdim * sdim)) / 2) + + fm = np.fft.fft2(mask) + pf = np.real(np.fft.ifft2(fm * np.fft.fft2(pr_high_res))) / np.sum(mask) + if np.sum(indmask) > 0: + pr_high_res[~indmask] = 1 + pf = pf / ( + np.real(np.fft.ifft2(fm * np.fft.fft2(pr_high_res))) + / np.sum(pr_high_res) + / len(fm) + ) + + pf[indmask] = np.nan + return pf + + else: + raise RuntimeError( + 'the smooth value does not exist, choose from this list {None,"Gaussian","tophat"}' + ) def _check_smooth_value(smooth): @@ -206,14 +212,14 @@ def downscale( elif type_smooth == 1: - P_agg = _smoothconv(P_u, orig_res) - r_agg = _smoothconv(r, orig_res) + P_agg = _balanced_spatial_average(P_u, orig_res=orig_res, kernel="Gaussian") + r_agg = _balanced_spatial_average(r, orig_res=orig_res, kernel="Gaussian") r *= P_agg / r_agg elif type_smooth == 2: - P_agg = _balanced_spatial_average(P_u, tophat) - r_agg = _balanced_spatial_average(r, tophat) + P_agg = _balanced_spatial_average(P_u, tophat=tophat, kernel="tophat") + r_agg = _balanced_spatial_average(r, tophat=tophat, kernel="tophat") r *= P_agg / r_agg if threshold is not None: From 388ba4560598153ec4125fc929c49a899e9d6c0e Mon Sep 17 00:00:00 2001 From: simonbeylat Date: Tue, 25 Oct 2022 14:32:31 +0200 Subject: [PATCH 04/21] add spectral merging from D'Onofrio et al. 2014 --- doc/source/references.bib | 11 +++++ pysteps/downscaling/rainfarm.py | 74 +++++++++++++++++++++++++++++---- 2 files changed, 78 insertions(+), 7 deletions(-) diff --git a/doc/source/references.bib b/doc/source/references.bib index 9f5eaf532..d4f638d94 100644 --- a/doc/source/references.bib +++ b/doc/source/references.bib @@ -366,4 +366,15 @@ @Article{Terzago2018 NUMBER = 11, PAGES = "2825--2840", DOI = "10.5194/nhess-18-2825-2018" +} + +@ARTICLE{DOnofrio2014, + TITLE = "Stochastic rainfall downscaling of climate models", + AUTHOR = "D'Onofrio, D and Palazzi, E and von Hardenberg, J and Provenzale, A and Calmanti, S", + JOURNAL = "J. Hydrometeorol.", + PUVLISHER = "American Meteorological Society", + VOLUME = 15, + NUMBER = 2, + PAGES = "830--843", + YEAR = 2014, } \ No newline at end of file diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index e6f83b0f3..1e54f2a50 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -22,6 +22,7 @@ import numpy as np from scipy.ndimage import convolve from scipy.ndimage import zoom +from pysteps.utils import spectral def _log_slope(log_k, log_power_spectrum): @@ -120,8 +121,48 @@ def _check_smooth_value(smooth): ) +def gaussianize(precip): + """ + Gaussianize field using rank ordering + + Parameters + ---------- + precip : matrix + matrix containing the input field to be Gaussianized. + + Returns + ------- + The Gaussianized field with the same dimensions as the input field. + + References + ---------- + :cite:`DOnofrio2014` + + """ + m = np.shape(precip)[0] + n = np.shape(precip)[1] + nn = m * n + ii = np.argsort(precip.reshape(nn)) + precip_gaussianize = np.zeros(nn) + precip_gaussianize[ii] = sorted(np.random.normal(0, 1, nn)) + precip_gaussianize = precip_gaussianize.reshape(m, n) + sd = np.std(precip_gaussianize) + if sd == 0: + sd = 1 + + precip_gaussianize /= sd + + return precip_gaussianize + + def downscale( - precip, ds_factor, alpha=None, threshold=None, return_alpha=False, smooth=None + precip, + ds_factor, + alpha=None, + threshold=None, + return_alpha=False, + smooth="Gaussian", + spectral_fusion=False, ): """ Downscale a rainfall field by increasing its spatial resolution by @@ -166,21 +207,25 @@ def downscale( """ type_smooth = _check_smooth_value(smooth) - orig_res = np.shape(precip)[1] - ki = np.fft.fftfreq(precip.shape[0]) - kj = np.fft.fftfreq(precip.shape[1]) + if spectral_fusion: + precip_bis = gaussianize(precip) + else: + precip_bis = precip + + ki = np.fft.fftfreq(precip_bis.shape[0]) + kj = np.fft.fftfreq(precip_bis.shape[1]) k_sqr = ki[:, None] ** 2 + kj[None, :] ** 2 k = np.sqrt(k_sqr) - ki_ds = np.fft.fftfreq(precip.shape[0] * ds_factor, d=1 / ds_factor) - kj_ds = np.fft.fftfreq(precip.shape[1] * ds_factor, d=1 / ds_factor) + ki_ds = np.fft.fftfreq(precip_bis.shape[0] * ds_factor, d=1 / ds_factor) + kj_ds = np.fft.fftfreq(precip_bis.shape[1] * ds_factor, d=1 / ds_factor) k_ds_sqr = ki_ds[:, None] ** 2 + kj_ds[None, :] ** 2 k_ds = np.sqrt(k_ds_sqr) if alpha is None: - fp = np.fft.fft2(precip) + fp = np.fft.fft2(precip_bis) fp_abs = abs(fp) log_power_spectrum = np.log(fp_abs**2) valid = (k != 0) & np.isfinite(log_power_spectrum) @@ -191,6 +236,21 @@ def downscale( warnings.simplefilter("ignore") fg *= np.sqrt(k_ds_sqr ** (-alpha / 2)) fg[0, 0] = 0 + + if spectral_fusion: + fp = np.fft.fft2(precip_bis) + kmax = int(precip_bis.shape[0] / 2) + nx = fg.shape[0] + Pk0 = spectral.rapsd(precip_bis)[kmax - 1] + gk0 = spectral.rapsd(np.fft.ifft2(fg).real / np.fft.ifft2(fg).real.std())[ + kmax - 1 + ] + fg *= Pk0 / gk0 + fg[0:kmax, 0:kmax] = fp[0:kmax, 0:kmax] + fg[nx - kmax :, 0:kmax] = fp[-kmax:, 0:kmax] + fg[0:kmax, nx - kmax :] = fp[0:kmax, -kmax:] + fg[nx - kmax :, nx - kmax :] = fp[-kmax:, -kmax:] + g = np.fft.ifft2(fg).real g /= g.std() r = np.exp(g) From 974f525c9cc3a59be92d48bc25d66cabbada548c Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Sun, 27 Nov 2022 15:46:30 +0100 Subject: [PATCH 05/21] Sort alphabetically --- doc/source/references.bib | 45 +++++++++++++++++++-------------------- 1 file changed, 22 insertions(+), 23 deletions(-) diff --git a/doc/source/references.bib b/doc/source/references.bib index d4f638d94..b78c58009 100644 --- a/doc/source/references.bib +++ b/doc/source/references.bib @@ -51,6 +51,17 @@ @ARTICLE{CRS2004 DOI = "10.1017/S1350482704001239" } +@ARTICLE{DOnofrio2014, + TITLE = "Stochastic rainfall downscaling of climate models", + AUTHOR = "D'Onofrio, D and Palazzi, E and von Hardenberg, J and Provenzale, A and Calmanti, S", + JOURNAL = "J. Hydrometeorol.", + PUVLISHER = "American Meteorological Society", + VOLUME = 15, + NUMBER = 2, + PAGES = "830--843", + YEAR = 2014, +} + @ARTICLE{EWWM2013, AUTHOR = "E. Ebert and L. Wilson and A. Weigel and M. Mittermaier and P. Nurmi and P. Gill and M. Göber and S. Joslyn and B. Brown and T. Fowler and A. Watkins", TITLE = "Progress and challenges in forecast verification", @@ -306,6 +317,17 @@ @ARTICLE{SPN2013 DOI = "10.1002/wrcr.20536" } +@Article{Terzago2018, + AUTHOR = "Terzago, S. and Palazzi, E. and von Hardenberg, J.", + TITLE = "Stochastic downscaling of precipitation in complex orography: a simple method to reproduce a realistic fine-scale climatology", + JOURNAL = "Natural Hazards and Earth System Sciences", + VOLUME = 18, + YEAR = 2018, + NUMBER = 11, + PAGES = "2825--2840", + DOI = "10.5194/nhess-18-2825-2018" +} + @ARTICLE{TRT2004, AUTHOR = "A. M. Hering and C. Morel and G. Galli and P. Ambrosetti and M. Boscacci", TITLE = "Nowcasting thunderstorms in the Alpine Region using a radar based adaptive thresholding scheme", @@ -355,26 +377,3 @@ @ARTICLE{ZR2009 YEAR = 2009, DOI = "10.1016/j.atmosres.2009.03.004" } - - -@Article{Terzago2018, -AUTHOR = "Terzago, S. and Palazzi, E. and von Hardenberg, J.", -TITLE = "Stochastic downscaling of precipitation in complex orography: a simple method to reproduce a realistic fine-scale climatology", -JOURNAL = "Natural Hazards and Earth System Sciences", -VOLUME = 18, -YEAR = 2018, -NUMBER = 11, -PAGES = "2825--2840", -DOI = "10.5194/nhess-18-2825-2018" -} - -@ARTICLE{DOnofrio2014, - TITLE = "Stochastic rainfall downscaling of climate models", - AUTHOR = "D'Onofrio, D and Palazzi, E and von Hardenberg, J and Provenzale, A and Calmanti, S", - JOURNAL = "J. Hydrometeorol.", - PUVLISHER = "American Meteorological Society", - VOLUME = 15, - NUMBER = 2, - PAGES = "830--843", - YEAR = 2014, -} \ No newline at end of file From 377764f54b22351e4824d9fbddab0c73907397ac Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Sun, 27 Nov 2022 16:05:47 +0100 Subject: [PATCH 06/21] Small fixes in the docstrings --- pysteps/downscaling/rainfarm.py | 38 +++++++++++++++++---------------- 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index 1e54f2a50..0107e7b2a 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -46,14 +46,14 @@ def _balanced_spatial_average(pr_high_res, tophat=None, orig_res=None, kernel="t """ Parameters ---------- - pr_high_res : matrix - matrix with the input field to smoothen, with dimensions ns*ns. - tophat : array_like - Array of weights - orig_res : int - original size of the precipitation field. - kernel : string - Choose the smoothing method + pr_high_res: array_like + Array containing the input precipitation field of shape (m, n). + tophat: array_like, optional + Array of weights, used only if kernel='tophat'. + orig_res: int, optional + Original size of the precipitation field, used only if kernel='Gaussian' + kernel : string, {None, "Gaussian", "tophat"} + The smoothing method. Returns ------- @@ -103,7 +103,7 @@ def _balanced_spatial_average(pr_high_res, tophat=None, orig_res=None, kernel="t else: raise RuntimeError( - 'the smooth value does not exist, choose from this list {None,"Gaussian","tophat"}' + 'the kernel value does not exist, choose from this list {None, "Gaussian", "tophat"}' ) @@ -117,22 +117,23 @@ def _check_smooth_value(smooth): return 2 else: raise RuntimeError( - 'the smooth value does not exist, choose from this list {None,"Gaussian","tophat"}' + 'the smooth value does not exist, choose from this list {None, "Gaussian", "tophat"}' ) def gaussianize(precip): """ - Gaussianize field using rank ordering + Gaussianize field using rank ordering. Parameters ---------- - precip : matrix - matrix containing the input field to be Gaussianized. + precip: array_like + Array containing the input field to be Gaussianized. Returns ------- - The Gaussianized field with the same dimensions as the input field. + precip_gaussianize: array_like + The Gaussianized field with the same shape as the input field. References ---------- @@ -170,7 +171,7 @@ def downscale( Parameters ---------- - precip: array-like + 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. @@ -184,11 +185,12 @@ def downscale( Set all values lower than the threshold to zero. return_alpha: bool, optional Whether to return the estimated spectral slope ``alpha``. - smooth: list, optional - Add the smoothing operatore from this list {None,"Gaussian","tophat"} + smooth: str, {None, "Gaussian", "tophat"} + The name of the smoothing operator. + Returns ------- - r: array-like + r: array_like Array of shape (m * ds_factor, n * ds_factor) containing the downscaled field. alpha: float From 694e9f5711d645b8a7a352fbd344247af331cd72 Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Sun, 27 Nov 2022 18:57:37 +0100 Subject: [PATCH 07/21] Refactor code --- pysteps/downscaling/rainfarm.py | 276 ++++++++++++-------------------- 1 file changed, 102 insertions(+), 174 deletions(-) diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index 0107e7b2a..0642bcbba 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -20,9 +20,31 @@ import warnings import numpy as np -from scipy.ndimage import convolve from scipy.ndimage import zoom -from pysteps.utils import spectral +from scipy.signal import convolve + + +def _gaussianize(precip): + """ + Gaussianize field using rank ordering as in :cite:`DOnofrio2014`. + """ + m, n = np.shape(precip) + nn = m * n + ii = np.argsort(precip.reshape(nn)) + precip_gaussianize = np.zeros(nn) + precip_gaussianize[ii] = sorted(np.random.normal(0, 1, nn)) + precip_gaussianize = precip_gaussianize.reshape(m, n) + sd = np.std(precip_gaussianize) + if sd == 0: + sd = 1 + return precip_gaussianize / sd + + +def _compute_freq_array(array, ds_factor=1): + freq_i = np.fft.fftfreq(array.shape[0] * ds_factor, d=1 / ds_factor) + freq_j = np.fft.fftfreq(array.shape[1] * ds_factor, d=1 / ds_factor) + freq_sqr = freq_i[:, None] ** 2 + freq_j[None, :] ** 2 + return np.sqrt(freq_sqr) def _log_slope(log_k, log_power_spectrum): @@ -31,129 +53,71 @@ def _log_slope(log_k, log_power_spectrum): lk_range = lk_max - lk_min lk_min += (1 / 6) * lk_range lk_max -= (1 / 6) * lk_range - selected = (lk_min <= log_k) & (log_k <= lk_max) lk_sel = log_k[selected] ps_sel = log_power_spectrum[selected] alpha = np.polyfit(lk_sel, ps_sel, 1)[0] alpha = -alpha - return alpha -def _balanced_spatial_average(pr_high_res, tophat=None, orig_res=None, kernel="tophat"): - - """ - Parameters - ---------- - pr_high_res: array_like - Array containing the input precipitation field of shape (m, n). - tophat: array_like, optional - Array of weights, used only if kernel='tophat'. - orig_res: int, optional - Original size of the precipitation field, used only if kernel='Gaussian' - kernel : string, {None, "Gaussian", "tophat"} - The smoothing method. - - Returns - ------- - The smoothened field. - - References - ---------- - :cite:`Terzago2018` +def _estimate_alpha(array, k): + fp = np.fft.fft2(array) + fp_abs = abs(fp) + log_power_spectrum = np.log(fp_abs**2) + valid = (k != 0) & np.isfinite(log_power_spectrum) + alpha = _log_slope(np.log(k[valid]), log_power_spectrum[valid]) + return alpha - """ - indmask = ~np.isfinite(pr_high_res) - pr_high_res[indmask] = 0.0 - - if kernel == "tophat": - ones = np.ones_like(pr_high_res) - return convolve(pr_high_res, tophat) / convolve(ones, tophat) - - elif kernel == "Gaussian": - ns = np.shape(pr_high_res)[1] - sdim = (ns / orig_res) / 2 - mask = np.zeros([ns, ns]) - for i in range(ns): - for j in range(ns): - - kx = i - ky = j - if i > (ns / 2): - kx = i - ns - if j > (ns / 2): - ky = j - ns - r2 = kx * kx + ky * ky - mask[i, j] = np.exp(-(r2 / (sdim * sdim)) / 2) - - fm = np.fft.fft2(mask) - pf = np.real(np.fft.ifft2(fm * np.fft.fft2(pr_high_res))) / np.sum(mask) - if np.sum(indmask) > 0: - pr_high_res[~indmask] = 1 - pf = pf / ( - np.real(np.fft.ifft2(fm * np.fft.fft2(pr_high_res))) - / np.sum(pr_high_res) - / len(fm) - ) - - pf[indmask] = np.nan - return pf +def _apply_spectral_fusion(array_low, array_high): + # TODO: implement same as in + # https://github.com/jhardenberg/RainFARM.jl/blob/master/src/rf/mergespec_spaceonly.jl + # https://github.com/jhardenberg/rainfarmr/blob/master/R/smoothconv.R + raise NotImplementedError - else: - raise RuntimeError( - 'the kernel value does not exist, choose from this list {None, "Gaussian", "tophat"}' - ) +def _compute_kernel_radius(ds_factor): + return int(round(ds_factor / np.sqrt(np.pi))) -def _check_smooth_value(smooth): - if smooth is None: - return 0 - elif smooth == "Gaussian": - return 1 - elif smooth == "tophat": - return 2 - else: - raise RuntimeError( - 'the smooth value does not exist, choose from this list {None, "Gaussian", "tophat"}' - ) +def _make_tophat_kernel(ds_factor): + """Compute 2d uniform (tophat) kernel""" + radius = _compute_kernel_radius(ds_factor) + (mx, my) = np.mgrid[-radius : radius + 0.01, -radius : radius + 0.01] + tophat = ((mx**2 + my**2) <= radius**2).astype(float) + return tophat / tophat.sum() -def gaussianize(precip): +def _make_gaussian_kernel(ds_factor): """ - Gaussianize field using rank ordering. - - Parameters - ---------- - precip: array_like - Array containing the input field to be Gaussianized. - - Returns - ------- - precip_gaussianize: array_like - The Gaussianized field with the same shape as the input field. + Compute 2d gaussian kernel + ref: https://github.com/scipy/scipy/blob/de80faf9d3480b9dbb9b888568b64499e0e70c19/scipy/ndimage/_filters.py#L179 + the smoothing sigma has width half a large pixel + """ + radius = _compute_kernel_radius(ds_factor) + sigma = ds_factor / 2 + sigma2 = sigma * sigma + x = np.arange(-radius, radius + 1) + kern1d = np.exp(-0.5 / sigma2 * x**2) + kern2d = np.outer(kern1d, kern1d) + return kern2d / kern2d.sum() - References - ---------- - :cite:`DOnofrio2014` - """ - m = np.shape(precip)[0] - n = np.shape(precip)[1] - nn = m * n - ii = np.argsort(precip.reshape(nn)) - precip_gaussianize = np.zeros(nn) - precip_gaussianize[ii] = sorted(np.random.normal(0, 1, nn)) - precip_gaussianize = precip_gaussianize.reshape(m, n) - sd = np.std(precip_gaussianize) - if sd == 0: - sd = 1 +def _balanced_spatial_average(array, kernel): + array = array.copy() + mask_valid = np.isfinite(array) + array[~mask_valid] = 0.0 + array_conv = convolve(array, kernel, mode="warp") + array_conv /= convolve(mask_valid, kernel, mode="warp") + array_conv[~mask_valid] = np.nan + return array_conv - precip_gaussianize /= sd - return precip_gaussianize +_make_kernel = dict() +_make_kernel["gaussian"] = _make_gaussian_kernel +_make_kernel["tophat"] = _make_tophat_kernel +_make_kernel["uniform"] = _make_tophat_kernel def downscale( @@ -162,12 +126,12 @@ def downscale( alpha=None, threshold=None, return_alpha=False, - smooth="Gaussian", + kernel_type="gaussian", spectral_fusion=False, ): """ - Downscale a rainfall field by increasing its spatial resolution by - a positive integer factor. + Downscale a rainfall field by increasing its spatial resolution by a positive + integer factor. Parameters ---------- @@ -185,12 +149,12 @@ def downscale( Set all values lower than the threshold to zero. return_alpha: bool, optional Whether to return the estimated spectral slope ``alpha``. - smooth: str, {None, "Gaussian", "tophat"} + kernel_type: {"gaussian", "uniform", "tophat"} The name of the smoothing operator. Returns ------- - r: array_like + precip_highres: ndarray Array of shape (m * ds_factor, n * ds_factor) containing the downscaled field. alpha: float @@ -208,86 +172,50 @@ def downscale( :cite:`Rebora2006` """ - type_smooth = _check_smooth_value(smooth) - orig_res = np.shape(precip)[1] if spectral_fusion: - precip_bis = gaussianize(precip) + precip_transformed = _gaussianize(precip) else: - precip_bis = precip + precip_transformed = precip - ki = np.fft.fftfreq(precip_bis.shape[0]) - kj = np.fft.fftfreq(precip_bis.shape[1]) - k_sqr = ki[:, None] ** 2 + kj[None, :] ** 2 - k = np.sqrt(k_sqr) - - ki_ds = np.fft.fftfreq(precip_bis.shape[0] * ds_factor, d=1 / ds_factor) - kj_ds = np.fft.fftfreq(precip_bis.shape[1] * ds_factor, d=1 / ds_factor) - k_ds_sqr = ki_ds[:, None] ** 2 + kj_ds[None, :] ** 2 - k_ds = np.sqrt(k_ds_sqr) + freq_array = _compute_freq_array(precip_transformed) + freq_array_highres = _compute_freq_array(precip_transformed, ds_factor) if alpha is None: - fp = np.fft.fft2(precip_bis) - fp_abs = abs(fp) - log_power_spectrum = np.log(fp_abs**2) - valid = (k != 0) & np.isfinite(log_power_spectrum) - alpha = _log_slope(np.log(k[valid]), log_power_spectrum[valid]) + alpha = _estimate_alpha(precip_transformed, freq_array) - fg = np.exp(complex(0, 1) * 2 * np.pi * np.random.rand(*k_ds.shape)) + white_noise_field = np.random.rand(*freq_array_highres.shape) + white_noise_field_complex = np.exp(complex(0, 1) * 2 * np.pi * white_noise_field) with warnings.catch_warnings(): warnings.simplefilter("ignore") - fg *= np.sqrt(k_ds_sqr ** (-alpha / 2)) - fg[0, 0] = 0 + noise_field_complex = white_noise_field_complex * np.sqrt( + freq_array_highres**-alpha + ) + noise_field_complex[0, 0] = 0 + noise_field = np.fft.ifft2(noise_field_complex).real + noise_field /= noise_field.std() + noise_field = np.exp(noise_field) if spectral_fusion: - fp = np.fft.fft2(precip_bis) - kmax = int(precip_bis.shape[0] / 2) - nx = fg.shape[0] - Pk0 = spectral.rapsd(precip_bis)[kmax - 1] - gk0 = spectral.rapsd(np.fft.ifft2(fg).real / np.fft.ifft2(fg).real.std())[ - kmax - 1 - ] - fg *= Pk0 / gk0 - fg[0:kmax, 0:kmax] = fp[0:kmax, 0:kmax] - fg[nx - kmax :, 0:kmax] = fp[-kmax:, 0:kmax] - fg[0:kmax, nx - kmax :] = fp[0:kmax, -kmax:] - fg[nx - kmax :, nx - kmax :] = fp[-kmax:, -kmax:] - - g = np.fft.ifft2(fg).real - g /= g.std() - r = np.exp(g) - - P_u = np.repeat(np.repeat(precip, ds_factor, axis=0), ds_factor, axis=1) - rad = int(round(ds_factor / np.sqrt(np.pi))) - (mx, my) = np.mgrid[-rad : rad + 0.01, -rad : rad + 0.01] - tophat = ((mx**2 + my**2) <= rad**2).astype(float) - tophat /= tophat.sum() - - if type_smooth == 0: - - P_agg = precip - r_agg = zoom(r, 1 / ds_factor, order=1) - factor = np.repeat( - np.repeat(P_agg / r_agg, ds_factor, axis=0), ds_factor, axis=1 + noise_field = _apply_spectral_fusion(precip_transformed, noise_field) + + try: + kernel = _make_kernel[kernel_type](ds_factor) + except KeyError: + raise ValueError( + f"kernel type '{kernel_type}' is invalid, " + f"available kernels: {list(_make_kernel)}" ) - r *= factor - - elif type_smooth == 1: - - P_agg = _balanced_spatial_average(P_u, orig_res=orig_res, kernel="Gaussian") - r_agg = _balanced_spatial_average(r, orig_res=orig_res, kernel="Gaussian") - r *= P_agg / r_agg - - elif type_smooth == 2: - P_agg = _balanced_spatial_average(P_u, tophat=tophat, kernel="tophat") - r_agg = _balanced_spatial_average(r, tophat=tophat, kernel="tophat") - r *= P_agg / r_agg + precip_zoomed = zoom(precip, ds_factor, order=1) + precip_smoothed = _balanced_spatial_average(precip_zoomed, kernel) + noise_smoothed = _balanced_spatial_average(noise_field, kernel) + precip_highres = noise_field * precip_smoothed / noise_smoothed if threshold is not None: - r[r < threshold] = 0 + precip_highres[precip_highres < threshold] = 0 if return_alpha: - return r, alpha + return precip_highres, alpha - return r + return precip_highres From 98db387dafecede0b01b6012325bf3a94de48416 Mon Sep 17 00:00:00 2001 From: simonbeylat Date: Tue, 6 Dec 2022 16:27:11 +0100 Subject: [PATCH 08/21] test smoothing option - ok, add first version of spectral fusion (still need to be tested and corrected) --- pysteps/downscaling/rainfarm.py | 50 +++++++++++++++++++++++++++++---- 1 file changed, 44 insertions(+), 6 deletions(-) diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index 0642bcbba..4e67ccf2e 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -71,10 +71,48 @@ def _estimate_alpha(array, k): def _apply_spectral_fusion(array_low, array_high): - # TODO: implement same as in - # https://github.com/jhardenberg/RainFARM.jl/blob/master/src/rf/mergespec_spaceonly.jl - # https://github.com/jhardenberg/rainfarmr/blob/master/R/smoothconv.R - raise NotImplementedError + + kmax = array_low.shape[0] // 2 + + (nax, nay) = np.shape(array_low) + (nx, ny) = np.shape(array_high) + + DFTr = np.fft.fft2(array_low) + DFTf = np.fft.fft2(array_high) + + DFTr = np.fft.fftshift(DFTr) + DFTf = np.fft.fftshift(DFTf) + + DFTfm = np.zeros([nx, nx]) + DFTr2 = np.zeros([nax + 1, nax + 1]) + + DFTr2[0:nax, 0:nax] = DFTr[:, :] + DFTr2[-1, 0:nax] = np.conj(DFTr[0, :]) + DFTr2[0:nax, -1] = np.conj(DFTr[:, 0]) + + kmax2 = kmax**2 + + ddx = 2 * np.pi / nax / 2 - 2 * np.pi / nx / 2 + + for i in range(nx): + for j in range(nx): + kx = -(nx // 2) + i - 1 + ky = -(nx // 2) + j - 1 + k2 = kx**2 + ky**2 + ir = (nax // 2) + 1 + kx + jr = (nax // 2) + 1 + ky + if k2 <= kmax2: + DFTfm[i, j] = DFTr2[ir - 1, jr - 1] * np.exp( + complex(0, ddx * kx - ddx * ky) + ) + else: + DFTfm[i, j] = DFTf[i, j] + + DFTfm = np.fft.ifftshift(DFTfm) + fm = np.fft.ifft2(DFTfm).real + fm /= fm.std() + fm = np.exp(fm) + return fm def _compute_kernel_radius(ds_factor): @@ -108,8 +146,8 @@ def _balanced_spatial_average(array, kernel): array = array.copy() mask_valid = np.isfinite(array) array[~mask_valid] = 0.0 - array_conv = convolve(array, kernel, mode="warp") - array_conv /= convolve(mask_valid, kernel, mode="warp") + array_conv = convolve(array, kernel, mode="same") + array_conv /= convolve(mask_valid, kernel, mode="same") array_conv[~mask_valid] = np.nan return array_conv From b8a8bebb5769863bb13324a2e46e29a4f59b3e1a Mon Sep 17 00:00:00 2001 From: blaiDoAr Date: Thu, 9 May 2024 11:55:37 +0200 Subject: [PATCH 09/21] Modify application of spectral fusion --- pysteps/downscaling/rainfarm.py | 117 +++++++++++++++++++++++--------- 1 file changed, 85 insertions(+), 32 deletions(-) diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index 4e67ccf2e..dea433692 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -68,51 +68,98 @@ def _estimate_alpha(array, k): valid = (k != 0) & np.isfinite(log_power_spectrum) alpha = _log_slope(np.log(k[valid]), log_power_spectrum[valid]) return alpha +#################################################### +####AFEGIM FUNCIO DE RADI MITJÀ DE POTÈNCIA ESPECTRAL: +################################################### +def fft2d_fromR(z): + ns = z.shape[0] + ny = z.shape[1] + if ns != ny: + raise ValueError("The input matrix must be square.") -def _apply_spectral_fusion(array_low, array_high): + ns2 = ns // 2 + nt = z.shape[2] if len(z.shape) > 2 else 1 + + if nt == 1: + z = np.expand_dims(z, axis=2) + + kx = np.tile(np.concatenate((np.arange(0, ns//2 + 1), np.arange(-ns//2 + 1, 0))), ns).reshape((ns, ns)) + km = np.sqrt(kx ** 2 + kx.T ** 2) + + fx = np.zeros(ns) + nn = np.zeros(ns) + + zf = np.abs(np.fft.fft2(z[:, :]) / (ns * ns)) ** 2 + zf0 = zf.copy() + zf[ns2, :] = zf0[ns2, :] / 2 + zf[:, ns2] = zf0[:, ns2] / 2 + for i in range(ns * ns): + ik = int(np.floor(km.flat[i] + 1.5)) + fx[ik] += zf.flat[i] + nn[ik] += 1 + + fx = fx[1:ns2 + 1] / nn[1:ns2 + 1] + + return fx +##################################### +##################################### + +def _apply_spectral_fusion(array_low, array_high,ds_factor): - kmax = array_low.shape[0] // 2 (nax, nay) = np.shape(array_low) (nx, ny) = np.shape(array_high) + nax2 = nax //2 + + ddx = np.pi*(1/nax - 1/nx) + + pstr = fft2d_fromR(array_high)*(nx*nx)**2 + pstra = fft2d_fromR(array_low)*(nax*nax)**2 + + c = pstra[nax2-1]/pstr[nax2-1] + array_high *= np.sqrt(c) + + DFTr = np.fft.fft2(array_low) DFTf = np.fft.fft2(array_high) - DFTr = np.fft.fftshift(DFTr) - DFTf = np.fft.fftshift(DFTf) - DFTfm = np.zeros([nx, nx]) - DFTr2 = np.zeros([nax + 1, nax + 1]) + DFTr2 = np.zeros((nx, nx), dtype=np.complex128) + DFTr2[0:nax2, 0:nax2] = DFTr[0:nax2, 0:nax2] + DFTr2[nx - nax2:nx, 0:nax2] = DFTr[nax2:2*nax2, 0:nax2] + DFTr2[0:nax2, nx - nax2:nx] = DFTr[0:nax2, nax2:2*nax2] + DFTr2[nx - nax2:nx, nx - nax2:nx] = DFTr[nax2:2*nax2, nax2:2*nax2] + + DFTr2[nax2, 0] = np.conj(DFTr2[nx - nax2, 0]) + DFTr2[0, nax2] = np.conj(DFTr2[0, nx - nax2]) - DFTr2[0:nax, 0:nax] = DFTr[:, :] - DFTr2[-1, 0:nax] = np.conj(DFTr[0, :]) - DFTr2[0:nax, -1] = np.conj(DFTr[:, 0]) - kmax2 = kmax**2 - ddx = 2 * np.pi / nax / 2 - 2 * np.pi / nx / 2 - for i in range(nx): - for j in range(nx): - kx = -(nx // 2) + i - 1 - ky = -(nx // 2) + j - 1 - k2 = kx**2 + ky**2 - ir = (nax // 2) + 1 + kx - jr = (nax // 2) + 1 + ky - if k2 <= kmax2: - DFTfm[i, j] = DFTr2[ir - 1, jr - 1] * np.exp( - complex(0, ddx * kx - ddx * ky) - ) - else: - DFTfm[i, j] = DFTf[i, j] + freq_array = _compute_freq_array(array_low) + freq_array_highres = _compute_freq_array(array_low, ds_factor) + + fax = np.max(freq_array) + + fx2 = freq_array_highres**2 + fax2 =fax**2 + + + fi = np.tile(np.fft.fftfreq(array_high.shape[0] , d=1 / ds_factor), nx).reshape((nx, nx)) + fj = fi.T + + DFTf = DFTf * (fx2 > fax2) + DFTr2 * (fx2 <= fax2) * np.exp(-1j * ddx * fi/fi[1,1] - 1j * ddx * fj/fj[1,1]) + + r = np.real(np.fft.ifftn(DFTf)) / len(DFTf) + + r /= r.std() + r = np.exp(r) + + return r + - DFTfm = np.fft.ifftshift(DFTfm) - fm = np.fft.ifft2(DFTfm).real - fm /= fm.std() - fm = np.exp(fm) - return fm def _compute_kernel_radius(ds_factor): @@ -231,11 +278,17 @@ def downscale( ) noise_field_complex[0, 0] = 0 noise_field = np.fft.ifft2(noise_field_complex).real - noise_field /= noise_field.std() - noise_field = np.exp(noise_field) if spectral_fusion: - noise_field = _apply_spectral_fusion(precip_transformed, noise_field) + + noise_field /= noise_field.shape[0] ** 2 + noise_field = np.exp(noise_field) + + noise_field = _apply_spectral_fusion(precip_transformed, noise_field, ds_factor) + + else: + noise_field /= noise_field.std() + noise_field = np.exp(noise_field) try: kernel = _make_kernel[kernel_type](ds_factor) From 4f959730a587fb3866928f34b76c7da180b83b55 Mon Sep 17 00:00:00 2001 From: jrmiro Date: Thu, 9 May 2024 13:54:01 +0200 Subject: [PATCH 10/21] Replace R radially averaged power spectrum function by rapsd PySTEPS function --- pysteps/downscaling/rainfarm.py | 54 +++++---------------------------- 1 file changed, 8 insertions(+), 46 deletions(-) diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index dea433692..41277a4d8 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -22,6 +22,7 @@ import numpy as np from scipy.ndimage import zoom from scipy.signal import convolve +from pysteps.utils.spectral import rapsd def _gaussianize(precip): @@ -68,59 +69,20 @@ def _estimate_alpha(array, k): valid = (k != 0) & np.isfinite(log_power_spectrum) alpha = _log_slope(np.log(k[valid]), log_power_spectrum[valid]) return alpha -#################################################### -####AFEGIM FUNCIO DE RADI MITJÀ DE POTÈNCIA ESPECTRAL: -################################################### -def fft2d_fromR(z): - ns = z.shape[0] - ny = z.shape[1] - if ns != ny: - raise ValueError("The input matrix must be square.") - ns2 = ns // 2 - nt = z.shape[2] if len(z.shape) > 2 else 1 - - if nt == 1: - z = np.expand_dims(z, axis=2) - - kx = np.tile(np.concatenate((np.arange(0, ns//2 + 1), np.arange(-ns//2 + 1, 0))), ns).reshape((ns, ns)) - km = np.sqrt(kx ** 2 + kx.T ** 2) - - fx = np.zeros(ns) - nn = np.zeros(ns) - - zf = np.abs(np.fft.fft2(z[:, :]) / (ns * ns)) ** 2 - zf0 = zf.copy() - zf[ns2, :] = zf0[ns2, :] / 2 - zf[:, ns2] = zf0[:, ns2] / 2 - for i in range(ns * ns): - ik = int(np.floor(km.flat[i] + 1.5)) - fx[ik] += zf.flat[i] - nn[ik] += 1 - - fx = fx[1:ns2 + 1] / nn[1:ns2 + 1] - - return fx -##################################### -##################################### - -def _apply_spectral_fusion(array_low, array_high,ds_factor): +def _apply_spectral_fusion(array_low, array_high, ds_factor): (nax, nay) = np.shape(array_low) (nx, ny) = np.shape(array_high) - nax2 = nax //2 + nax2 = nax //2 - ddx = np.pi*(1/nax - 1/nx) - - pstr = fft2d_fromR(array_high)*(nx*nx)**2 - pstra = fft2d_fromR(array_low)*(nax*nax)**2 - - c = pstra[nax2-1]/pstr[nax2-1] - array_high *= np.sqrt(c) + array_low_k0 = rapsd(array_low, fft_method=np.fft)[nax2 - 1] * nax ** 2 + array_high_k0 = rapsd(array_high, fft_method=np.fft)[nax2 - 1] * nx ** 2 + array_high *= np.sqrt(array_low_k0 / array_high_k0) DFTr = np.fft.fft2(array_low) DFTf = np.fft.fft2(array_high) @@ -150,6 +112,8 @@ def _apply_spectral_fusion(array_low, array_high,ds_factor): fi = np.tile(np.fft.fftfreq(array_high.shape[0] , d=1 / ds_factor), nx).reshape((nx, nx)) fj = fi.T + ddx = np.pi*(1/nax - 1/nx) + DFTf = DFTf * (fx2 > fax2) + DFTr2 * (fx2 <= fax2) * np.exp(-1j * ddx * fi/fi[1,1] - 1j * ddx * fj/fj[1,1]) r = np.real(np.fft.ifftn(DFTf)) / len(DFTf) @@ -160,8 +124,6 @@ def _apply_spectral_fusion(array_low, array_high,ds_factor): return r - - def _compute_kernel_radius(ds_factor): return int(round(ds_factor / np.sqrt(np.pi))) From f838b649c4fd22ab68dd311d11ff3c0f01587fe6 Mon Sep 17 00:00:00 2001 From: ecasellas Date: Fri, 10 May 2024 08:51:37 +0200 Subject: [PATCH 11/21] Include freq parameters into _apply_spectral_fusion and change variable names --- pysteps/downscaling/rainfarm.py | 66 ++++++++++++++------------------- 1 file changed, 27 insertions(+), 39 deletions(-) diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index 41277a4d8..fee7e647d 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -71,57 +71,45 @@ def _estimate_alpha(array, k): return alpha -def _apply_spectral_fusion(array_low, array_high, ds_factor): +def _apply_spectral_fusion(array_low, array_high, freq_array_low, freq_array_high, ds_factor): + nax, _ = np.shape(array_low) + nx, _ = np.shape(array_high) - (nax, nay) = np.shape(array_low) - (nx, ny) = np.shape(array_high) - - nax2 = nax //2 + nax2 = nax // 2 array_low_k0 = rapsd(array_low, fft_method=np.fft)[nax2 - 1] * nax ** 2 array_high_k0 = rapsd(array_high, fft_method=np.fft)[nax2 - 1] * nx ** 2 array_high *= np.sqrt(array_low_k0 / array_high_k0) - DFTr = np.fft.fft2(array_low) - DFTf = np.fft.fft2(array_high) - - - DFTr2 = np.zeros((nx, nx), dtype=np.complex128) - DFTr2[0:nax2, 0:nax2] = DFTr[0:nax2, 0:nax2] - DFTr2[nx - nax2:nx, 0:nax2] = DFTr[nax2:2*nax2, 0:nax2] - DFTr2[0:nax2, nx - nax2:nx] = DFTr[0:nax2, nax2:2*nax2] - DFTr2[nx - nax2:nx, nx - nax2:nx] = DFTr[nax2:2*nax2, nax2:2*nax2] - - DFTr2[nax2, 0] = np.conj(DFTr2[nx - nax2, 0]) - DFTr2[0, nax2] = np.conj(DFTr2[0, nx - nax2]) - - + fft_array_low = np.fft.fft2(array_low) + fft_array_high = np.fft.fft2(array_high) + fft_merged = np.zeros((nx, nx), dtype=np.complex128) + fft_merged[0:nax2, 0:nax2] = fft_array_low[0:nax2, 0:nax2] + fft_merged[nx - nax2:nx, 0:nax2] = fft_array_low[nax2:2*nax2, 0:nax2] + fft_merged[0:nax2, nx - nax2:nx] = fft_array_low[0:nax2, nax2:2*nax2] + fft_merged[nx - nax2:nx, nx - nax2:nx] = fft_array_low[nax2:2*nax2, nax2:2*nax2] - freq_array = _compute_freq_array(array_low) - freq_array_highres = _compute_freq_array(array_low, ds_factor) + fft_merged[nax2, 0] = np.conj(fft_merged[nx - nax2, 0]) + fft_merged[0, nax2] = np.conj(fft_merged[0, nx - nax2]) - fax = np.max(freq_array) + freq_i = np.tile(np.fft.fftfreq(array_high.shape[0], d=1 / ds_factor), nx).reshape((nx, nx)) + freq_j = freq_i.T - fx2 = freq_array_highres**2 - fax2 =fax**2 + ddx = np.pi * (1 / nax - 1 / nx) / np.abs(freq_i[0, 1] - freq_i[0, 0]) + fx2 = freq_array_high ** 2 + fax2 = np.max(freq_array_low) ** 2 + fft_merged = fft_array_high * (fx2 > fax2) + fft_merged * (fx2 <= fax2) * np.exp(-1j * ddx * freq_i - 1j * ddx * freq_j) - fi = np.tile(np.fft.fftfreq(array_high.shape[0] , d=1 / ds_factor), nx).reshape((nx, nx)) - fj = fi.T + merged = np.real(np.fft.ifftn(fft_merged)) / len(fft_merged) - ddx = np.pi*(1/nax - 1/nx) + merged /= merged.std() + merged = np.exp(merged) - DFTf = DFTf * (fx2 > fax2) + DFTr2 * (fx2 <= fax2) * np.exp(-1j * ddx * fi/fi[1,1] - 1j * ddx * fj/fj[1,1]) - - r = np.real(np.fft.ifftn(DFTf)) / len(DFTf) - - r /= r.std() - r = np.exp(r) - - return r + return merged def _compute_kernel_radius(ds_factor): @@ -236,18 +224,18 @@ def downscale( with warnings.catch_warnings(): warnings.simplefilter("ignore") noise_field_complex = white_noise_field_complex * np.sqrt( - freq_array_highres**-alpha + freq_array_highres**-alpha / 2 ) noise_field_complex[0, 0] = 0 noise_field = np.fft.ifft2(noise_field_complex).real if spectral_fusion: - noise_field /= noise_field.shape[0] ** 2 noise_field = np.exp(noise_field) - noise_field = _apply_spectral_fusion(precip_transformed, noise_field, ds_factor) - + noise_field = _apply_spectral_fusion( + precip_transformed, noise_field, freq_array, freq_array_highres, ds_factor + ) else: noise_field /= noise_field.std() noise_field = np.exp(noise_field) From 3ce2d77dc8c4d19b0b91ad63816a3d061baba366 Mon Sep 17 00:00:00 2001 From: ecasellas Date: Fri, 10 May 2024 11:50:22 +0200 Subject: [PATCH 12/21] Adapt downscale function to spectral_fusion --- pysteps/downscaling/rainfarm.py | 55 ++++++++++++++++++++------------- 1 file changed, 34 insertions(+), 21 deletions(-) diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index fee7e647d..287d928a3 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -23,6 +23,7 @@ from scipy.ndimage import zoom from scipy.signal import convolve from pysteps.utils.spectral import rapsd +from pysteps.utils.dimension import aggregate_fields def _gaussianize(precip): @@ -71,15 +72,17 @@ def _estimate_alpha(array, k): return alpha -def _apply_spectral_fusion(array_low, array_high, freq_array_low, freq_array_high, ds_factor): +def _apply_spectral_fusion( + array_low, array_high, freq_array_low, freq_array_high, ds_factor +): nax, _ = np.shape(array_low) nx, _ = np.shape(array_high) - nax2 = nax // 2 + k0 = nax // 2 - array_low_k0 = rapsd(array_low, fft_method=np.fft)[nax2 - 1] * nax ** 2 - array_high_k0 = rapsd(array_high, fft_method=np.fft)[nax2 - 1] * nx ** 2 + array_low_k0 = rapsd(array_low, fft_method=np.fft)[k0 - 1] * nax**2 + array_high_k0 = rapsd(array_high, fft_method=np.fft)[k0 - 1] * nx**2 array_high *= np.sqrt(array_low_k0 / array_high_k0) @@ -87,22 +90,26 @@ def _apply_spectral_fusion(array_low, array_high, freq_array_low, freq_array_hig fft_array_high = np.fft.fft2(array_high) fft_merged = np.zeros((nx, nx), dtype=np.complex128) - fft_merged[0:nax2, 0:nax2] = fft_array_low[0:nax2, 0:nax2] - fft_merged[nx - nax2:nx, 0:nax2] = fft_array_low[nax2:2*nax2, 0:nax2] - fft_merged[0:nax2, nx - nax2:nx] = fft_array_low[0:nax2, nax2:2*nax2] - fft_merged[nx - nax2:nx, nx - nax2:nx] = fft_array_low[nax2:2*nax2, nax2:2*nax2] + fft_merged[0:k0, 0:k0] = fft_array_low[0:k0, 0:k0] + fft_merged[nx - k0:nx, 0:k0] = fft_array_low[k0:2 * k0, 0:k0] + fft_merged[0:k0, nx - k0:nx] = fft_array_low[0:k0, k0:2 * k0] + fft_merged[nx - k0:nx, nx - k0:nx] = fft_array_low[k0:2 * k0, k0:2 * k0] - fft_merged[nax2, 0] = np.conj(fft_merged[nx - nax2, 0]) - fft_merged[0, nax2] = np.conj(fft_merged[0, nx - nax2]) + fft_merged[k0, 0] = np.conj(fft_merged[nx - k0, 0]) + fft_merged[0, k0] = np.conj(fft_merged[0, nx - k0]) - freq_i = np.tile(np.fft.fftfreq(array_high.shape[0], d=1 / ds_factor), nx).reshape((nx, nx)) + freq_i = np.tile(np.fft.fftfreq(array_high.shape[0], d=1 / ds_factor), nx).reshape( + (nx, nx) + ) freq_j = freq_i.T ddx = np.pi * (1 / nax - 1 / nx) / np.abs(freq_i[0, 1] - freq_i[0, 0]) - fx2 = freq_array_high ** 2 - fax2 = np.max(freq_array_low) ** 2 + fx2 = freq_array_high**2 + fax2 = freq_array_low[k0, k0] ** 2 - fft_merged = fft_array_high * (fx2 > fax2) + fft_merged * (fx2 <= fax2) * np.exp(-1j * ddx * freq_i - 1j * ddx * freq_j) + fft_merged = fft_array_high * (fx2 > fax2) + fft_merged * (fx2 <= fax2) * np.exp( + -1j * ddx * freq_i - 1j * ddx * freq_j + ) merged = np.real(np.fft.ifftn(fft_merged)) / len(fft_merged) @@ -119,7 +126,7 @@ def _compute_kernel_radius(ds_factor): def _make_tophat_kernel(ds_factor): """Compute 2d uniform (tophat) kernel""" radius = _compute_kernel_radius(ds_factor) - (mx, my) = np.mgrid[-radius : radius + 0.01, -radius : radius + 0.01] + (mx, my) = np.mgrid[-radius:radius + 0.01, -radius:radius + 0.01] tophat = ((mx**2 + my**2) <= radius**2).astype(float) return tophat / tophat.sum() @@ -235,8 +242,8 @@ def downscale( noise_field = _apply_spectral_fusion( precip_transformed, noise_field, freq_array, freq_array_highres, ds_factor - ) - else: + ) + else: noise_field /= noise_field.std() noise_field = np.exp(noise_field) @@ -248,10 +255,16 @@ def downscale( f"available kernels: {list(_make_kernel)}" ) - precip_zoomed = zoom(precip, ds_factor, order=1) - precip_smoothed = _balanced_spatial_average(precip_zoomed, kernel) - noise_smoothed = _balanced_spatial_average(noise_field, kernel) - precip_highres = noise_field * precip_smoothed / noise_smoothed + noise_lowres = aggregate_fields(noise_field, ds_factor, axis=(0, 1)) + + ca = precip / noise_lowres + + cai = np.repeat(np.repeat(ca, ds_factor, axis=0), ds_factor, axis=1) + + precip_highres = noise_field * cai + + if kernel is not None: + precip_highres = _balanced_spatial_average(precip_highres, kernel) if threshold is not None: precip_highres[precip_highres < threshold] = 0 From be1bb62b1d877b4af8a12a610200843ae03335c8 Mon Sep 17 00:00:00 2001 From: ecasellas Date: Fri, 10 May 2024 13:19:36 +0200 Subject: [PATCH 13/21] Correct smoothing application in downscale function --- pysteps/downscaling/rainfarm.py | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index 287d928a3..4472d52f5 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -247,24 +247,26 @@ def downscale( noise_field /= noise_field.std() noise_field = np.exp(noise_field) - try: - kernel = _make_kernel[kernel_type](ds_factor) - except KeyError: - raise ValueError( - f"kernel type '{kernel_type}' is invalid, " - f"available kernels: {list(_make_kernel)}" - ) - noise_lowres = aggregate_fields(noise_field, ds_factor, axis=(0, 1)) - ca = precip / noise_lowres + precip = np.kron(precip, np.ones((ds_factor, ds_factor))) + noise_lowres = np.kron(noise_lowres, np.ones((ds_factor, ds_factor))) + + if kernel_type is not None: + try: + kernel = _make_kernel[kernel_type](ds_factor) + except KeyError: + raise ValueError( + f"kernel type '{kernel_type}' is invalid, " + f"available kernels: {list(_make_kernel)}" + ) - cai = np.repeat(np.repeat(ca, ds_factor, axis=0), ds_factor, axis=1) + precip = _balanced_spatial_average(precip, kernel) + noise_lowres = _balanced_spatial_average(noise_lowres, kernel) - precip_highres = noise_field * cai + norm_k0 = precip / noise_lowres - if kernel is not None: - precip_highres = _balanced_spatial_average(precip_highres, kernel) + precip_highres = noise_field * norm_k0 if threshold is not None: precip_highres[precip_highres < threshold] = 0 From d500a06215483a0fb7fc5d120d2ceb3ca9965df8 Mon Sep 17 00:00:00 2001 From: ecasellas Date: Tue, 14 May 2024 07:48:46 +0200 Subject: [PATCH 14/21] Remove unused zoom and reduce code --- pysteps/downscaling/rainfarm.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index 4472d52f5..98aceedc5 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -20,7 +20,6 @@ import warnings import numpy as np -from scipy.ndimage import zoom from scipy.signal import convolve from pysteps.utils.spectral import rapsd from pysteps.utils.dimension import aggregate_fields @@ -113,9 +112,6 @@ def _apply_spectral_fusion( merged = np.real(np.fft.ifftn(fft_merged)) / len(fft_merged) - merged /= merged.std() - merged = np.exp(merged) - return merged @@ -165,7 +161,7 @@ def _balanced_spatial_average(array, kernel): def downscale( precip, ds_factor, - alpha=None, + alpha: None, threshold=None, return_alpha=False, kernel_type="gaussian", @@ -243,9 +239,9 @@ def downscale( noise_field = _apply_spectral_fusion( precip_transformed, noise_field, freq_array, freq_array_highres, ds_factor ) - else: - noise_field /= noise_field.std() - noise_field = np.exp(noise_field) + + noise_field /= noise_field.std() + noise_field = np.exp(noise_field) noise_lowres = aggregate_fields(noise_field, ds_factor, axis=(0, 1)) From bc4943ee1d87f9f866f4a5202b8bc62512a511e3 Mon Sep 17 00:00:00 2001 From: blaiDoAr Date: Tue, 14 May 2024 12:26:08 +0200 Subject: [PATCH 15/21] Add aggregation test and adapt docstring --- pysteps/downscaling/rainfarm.py | 13 ++-- pysteps/tests/test_downscaling_rainfarm.py | 87 ++++++++++++++++++++-- 2 files changed, 88 insertions(+), 12 deletions(-) diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index 98aceedc5..de666c0e3 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -4,7 +4,7 @@ ============================ Implementation of the RainFARM stochastic downscaling method as described in -:cite:`Rebora2006`. +:cite:`Rebora2006` and :cite:`DOnofrio2014`. 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 @@ -164,7 +164,7 @@ def downscale( alpha: None, threshold=None, return_alpha=False, - kernel_type="gaussian", + kernel_type=None, spectral_fusion=False, ): """ @@ -187,8 +187,10 @@ def downscale( Set all values lower than the threshold to zero. return_alpha: bool, optional Whether to return the estimated spectral slope ``alpha``. - kernel_type: {"gaussian", "uniform", "tophat"} - The name of the smoothing operator. + kernel_type: {None, "gaussian", "uniform", "tophat"} + The name of the smoothing operator. If None no smoothing is applied. + spectral_fusion: bool, optional + Whether to apply spectral merging as in :cite:`DOnofrio2014`. Returns ------- @@ -203,11 +205,12 @@ def downscale( Currently, the pysteps implementation of RainFARM only covers spatial downscaling. That is, it can improve the spatial resolution of a rainfall field. However, unlike the original algorithm from Rebora et al. (2006), it cannot downscale the temporal - dimension. + dimension. It implements spectral merging from D'Onofrio et al. (2014). References ---------- :cite:`Rebora2006` + :cite:`DOnofrio2014` """ diff --git a/pysteps/tests/test_downscaling_rainfarm.py b/pysteps/tests/test_downscaling_rainfarm.py index ec2a3fbb8..d1f42a6da 100644 --- a/pysteps/tests/test_downscaling_rainfarm.py +++ b/pysteps/tests/test_downscaling_rainfarm.py @@ -1,10 +1,10 @@ # -*- coding: utf-8 -*- import pytest - +import numpy as np from pysteps import downscaling from pysteps.tests.helpers import get_precipitation_fields -from pysteps.utils import aggregate_fields_space, square_domain +from pysteps.utils import aggregate_fields_space, square_domain, aggregate_fields @pytest.fixture(scope="module") @@ -17,12 +17,32 @@ def data(): 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)] +rainfarm_arg_names = ( + "alpha", + "ds_factor", + "threshold", + "return_alpha", + "spectral_fusion", + "kernel_type", +) +rainfarm_arg_values = [ + (1.0, 1, 0, False, False, None), + (1, 2, 0, False, False, "gaussian"), + (1, 4, 0, False, False, "tophat"), + (1, 4, 0, False, True, "uniform"), +] @pytest.mark.parametrize(rainfarm_arg_names, rainfarm_arg_values) -def test_rainfarm_shape(data, alpha, ds_factor, threshold, return_alpha): +def test_rainfarm_shape( + data, + alpha, + ds_factor, + threshold, + return_alpha, + spectral_fusion, + kernel_type, +): """Test that the output of rainfarm is consistent with the downscaling factor.""" precip, metadata = data window = metadata["xpixelsize"] * ds_factor @@ -35,6 +55,8 @@ def test_rainfarm_shape(data, alpha, ds_factor, threshold, return_alpha): ds_factor=ds_factor, threshold=threshold, return_alpha=return_alpha, + spectral_fusion=spectral_fusion, + kernel_type=kernel_type, ) assert precip_hr.ndim == precip.ndim @@ -42,11 +64,60 @@ def test_rainfarm_shape(data, alpha, ds_factor, threshold, return_alpha): assert precip_hr.shape[1] == precip.shape[1] -rainfarm_arg_values = [(1.0, 2, 0, True), (None, 2, 0, True)] +rainfarm_arg_values = [ + (1.0, 1, 0, False, False, None), + (1, 2, 0, False, False, None), + (1, 4, 0, False, False, None), + (1, 4, 0, False, True, None), +] + + +@pytest.mark.parametrize(rainfarm_arg_names, rainfarm_arg_values) +def test_rainfarm_aggregate( + data, + alpha, + ds_factor, + threshold, + return_alpha, + spectral_fusion, + kernel_type, +): + """Test that the output of rainfarm is equal to original when aggregated.""" + 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=alpha, + ds_factor=ds_factor, + threshold=threshold, + return_alpha=return_alpha, + spectral_fusion=spectral_fusion, + kernel_type=kernel_type, + ) + precip_low = aggregate_fields(precip_hr, ds_factor, axis=(0, 1)) + precip_lr[precip_lr < threshold] = 0.0 + + print(np.max(precip_low), np.max(precip_lr)) + + np.testing.assert_array_almost_equal(precip_lr, precip_low) + + +rainfarm_arg_values = [(1.0, 2, 0, True, False, None), (None, 2, 0, True, True, None)] @pytest.mark.parametrize(rainfarm_arg_names, rainfarm_arg_values) -def test_rainfarm_alpha(data, alpha, ds_factor, threshold, return_alpha): +def test_rainfarm_alpha( + data, + alpha, + ds_factor, + threshold, + return_alpha, + spectral_fusion, + kernel_type, +): """Test that rainfarm computes and returns alpha.""" precip, metadata = data window = metadata["xpixelsize"] * ds_factor @@ -59,6 +130,8 @@ def test_rainfarm_alpha(data, alpha, ds_factor, threshold, return_alpha): ds_factor=ds_factor, threshold=threshold, return_alpha=return_alpha, + spectral_fusion=spectral_fusion, + kernel_type=kernel_type, ) assert len(precip_hr) == 2 From 08cff0d120e36e4bdcac4c225966d328082f5027 Mon Sep 17 00:00:00 2001 From: blaiDoAr Date: Tue, 14 May 2024 12:34:00 +0200 Subject: [PATCH 16/21] Remove print --- pysteps/tests/test_downscaling_rainfarm.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/pysteps/tests/test_downscaling_rainfarm.py b/pysteps/tests/test_downscaling_rainfarm.py index d1f42a6da..884270d09 100644 --- a/pysteps/tests/test_downscaling_rainfarm.py +++ b/pysteps/tests/test_downscaling_rainfarm.py @@ -100,8 +100,6 @@ def test_rainfarm_aggregate( precip_low = aggregate_fields(precip_hr, ds_factor, axis=(0, 1)) precip_lr[precip_lr < threshold] = 0.0 - print(np.max(precip_low), np.max(precip_lr)) - np.testing.assert_array_almost_equal(precip_lr, precip_low) From 2974eaf8e19c94bbdaabb94ccc4432e10466eaa8 Mon Sep 17 00:00:00 2001 From: ecasellas Date: Tue, 14 May 2024 12:56:59 +0200 Subject: [PATCH 17/21] Set None as default alpha value --- pysteps/downscaling/rainfarm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index de666c0e3..8a2a4ac20 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -161,7 +161,7 @@ def _balanced_spatial_average(array, kernel): def downscale( precip, ds_factor, - alpha: None, + alpha=None, threshold=None, return_alpha=False, kernel_type=None, From f6f1b4c4fe34a36e0b804e01f30f38802d4c6fe5 Mon Sep 17 00:00:00 2001 From: ecasellas Date: Fri, 24 May 2024 13:27:05 +0200 Subject: [PATCH 18/21] Correct alpha value synthetic field amplitude --- pysteps/downscaling/rainfarm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index 8a2a4ac20..1d1033f18 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -230,7 +230,7 @@ def downscale( with warnings.catch_warnings(): warnings.simplefilter("ignore") noise_field_complex = white_noise_field_complex * np.sqrt( - freq_array_highres**-alpha / 2 + freq_array_highres**-alpha ) noise_field_complex[0, 0] = 0 noise_field = np.fft.ifft2(noise_field_complex).real From 51561706b83c7c0b0c44c6ca7b70121f1ea564ec Mon Sep 17 00:00:00 2001 From: ecasellas Date: Mon, 27 May 2024 07:50:23 +0200 Subject: [PATCH 19/21] Run black formatter on rainfarm.py --- pysteps/downscaling/rainfarm.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index 1d1033f18..69526eb50 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -90,9 +90,9 @@ def _apply_spectral_fusion( fft_merged = np.zeros((nx, nx), dtype=np.complex128) fft_merged[0:k0, 0:k0] = fft_array_low[0:k0, 0:k0] - fft_merged[nx - k0:nx, 0:k0] = fft_array_low[k0:2 * k0, 0:k0] - fft_merged[0:k0, nx - k0:nx] = fft_array_low[0:k0, k0:2 * k0] - fft_merged[nx - k0:nx, nx - k0:nx] = fft_array_low[k0:2 * k0, k0:2 * k0] + fft_merged[nx - k0 : nx, 0:k0] = fft_array_low[k0 : 2 * k0, 0:k0] + fft_merged[0:k0, nx - k0 : nx] = fft_array_low[0:k0, k0 : 2 * k0] + fft_merged[nx - k0 : nx, nx - k0 : nx] = fft_array_low[k0 : 2 * k0, k0 : 2 * k0] fft_merged[k0, 0] = np.conj(fft_merged[nx - k0, 0]) fft_merged[0, k0] = np.conj(fft_merged[0, nx - k0]) @@ -122,7 +122,7 @@ def _compute_kernel_radius(ds_factor): def _make_tophat_kernel(ds_factor): """Compute 2d uniform (tophat) kernel""" radius = _compute_kernel_radius(ds_factor) - (mx, my) = np.mgrid[-radius:radius + 0.01, -radius:radius + 0.01] + (mx, my) = np.mgrid[-radius : radius + 0.01, -radius : radius + 0.01] tophat = ((mx**2 + my**2) <= radius**2).astype(float) return tophat / tophat.sum() @@ -189,7 +189,7 @@ def downscale( Whether to return the estimated spectral slope ``alpha``. kernel_type: {None, "gaussian", "uniform", "tophat"} The name of the smoothing operator. If None no smoothing is applied. - spectral_fusion: bool, optional + spectral_fusion: bool, optional Whether to apply spectral merging as in :cite:`DOnofrio2014`. Returns @@ -205,7 +205,7 @@ def downscale( Currently, the pysteps implementation of RainFARM only covers spatial downscaling. That is, it can improve the spatial resolution of a rainfall field. However, unlike the original algorithm from Rebora et al. (2006), it cannot downscale the temporal - dimension. It implements spectral merging from D'Onofrio et al. (2014). + dimension. It implements spectral merging from D'Onofrio et al. (2014). References ---------- From c38f4e11b6801744ac97486ef31f586465c3c91b Mon Sep 17 00:00:00 2001 From: Daniele Nerini Date: Sun, 9 Jun 2024 11:55:53 +0200 Subject: [PATCH 20/21] Some refactoring to improve readability --- pysteps/downscaling/rainfarm.py | 141 ++++++++++++++++++++++---------- 1 file changed, 97 insertions(+), 44 deletions(-) diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index 69526eb50..d9906390d 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -42,6 +42,9 @@ def _gaussianize(precip): def _compute_freq_array(array, ds_factor=1): + """ + Compute the frequency array following a given downscaling factor. + """ freq_i = np.fft.fftfreq(array.shape[0] * ds_factor, d=1 / ds_factor) freq_j = np.fft.fftfreq(array.shape[1] * ds_factor, d=1 / ds_factor) freq_sqr = freq_i[:, None] ** 2 + freq_j[None, :] ** 2 @@ -49,6 +52,10 @@ def _compute_freq_array(array, ds_factor=1): def _log_slope(log_k, log_power_spectrum): + """ + Calculate the log-slope of the power spectrum given an array of logarithmic wavenumbers + and an array of logarithmic power spectrum values. + """ lk_min = log_k.min() lk_max = log_k.max() lk_range = lk_max - lk_min @@ -63,6 +70,9 @@ def _log_slope(log_k, log_power_spectrum): def _estimate_alpha(array, k): + """ + Estimate the alpha parameter using the power spectrum of the input array. + """ fp = np.fft.fft2(array) fp_abs = abs(fp) log_power_spectrum = np.log(fp_abs**2) @@ -71,45 +81,82 @@ def _estimate_alpha(array, k): return alpha +def _compute_noise_field(freq_array_highres, alpha): + """ + Compute a field of correlated noise field using the given frequency array and alpha + value. + """ + white_noise_field = np.random.rand(*freq_array_highres.shape) + white_noise_field_complex = np.exp(complex(0, 1) * 2 * np.pi * white_noise_field) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + noise_field_complex = white_noise_field_complex * np.sqrt( + freq_array_highres**-alpha + ) + noise_field_complex[0, 0] = 0 + return np.fft.ifft2(noise_field_complex).real + + def _apply_spectral_fusion( array_low, array_high, freq_array_low, freq_array_high, ds_factor ): + """ + Apply spectral fusion to merge two arrays in the frequency domain. + """ + + # Validate inputs + if array_low.shape != freq_array_low.shape: + raise ValueError("Shape of array_low must match shape of freq_array_low.") + if array_high.shape != freq_array_high.shape: + raise ValueError("Shape of array_high must match shape of freq_array_high.") nax, _ = np.shape(array_low) nx, _ = np.shape(array_high) - k0 = nax // 2 - array_low_k0 = rapsd(array_low, fft_method=np.fft)[k0 - 1] * nax**2 - array_high_k0 = rapsd(array_high, fft_method=np.fft)[k0 - 1] * nx**2 + # Calculate power spectral density at specific frequency + def compute_psd(array, fft_size): + return rapsd(array, fft_method=np.fft)[k0 - 1] * fft_size**2 + + psd_low = compute_psd(array_low, nax) + psd_high = compute_psd(array_high, nx) - array_high *= np.sqrt(array_low_k0 / array_high_k0) + # Normalize high-resolution array + normalization_factor = np.sqrt(psd_low / psd_high) + array_high *= normalization_factor - fft_array_low = np.fft.fft2(array_low) - fft_array_high = np.fft.fft2(array_high) + # Perform FFT on both arrays + fft_low = np.fft.fft2(array_low) + fft_high = np.fft.fft2(array_high) - fft_merged = np.zeros((nx, nx), dtype=np.complex128) - fft_merged[0:k0, 0:k0] = fft_array_low[0:k0, 0:k0] - fft_merged[nx - k0 : nx, 0:k0] = fft_array_low[k0 : 2 * k0, 0:k0] - fft_merged[0:k0, nx - k0 : nx] = fft_array_low[0:k0, k0 : 2 * k0] - fft_merged[nx - k0 : nx, nx - k0 : nx] = fft_array_low[k0 : 2 * k0, k0 : 2 * k0] + # Initialize the merged FFT array with low-resolution data + fft_merged = np.zeros_like(fft_high, dtype=np.complex128) + fft_merged[0:k0, 0:k0] = fft_low[0:k0, 0:k0] + fft_merged[nx - k0 : nx, 0:k0] = fft_low[k0 : 2 * k0, 0:k0] + fft_merged[0:k0, nx - k0 : nx] = fft_low[0:k0, k0 : 2 * k0] + fft_merged[nx - k0 : nx, nx - k0 : nx] = fft_low[k0 : 2 * k0, k0 : 2 * k0] fft_merged[k0, 0] = np.conj(fft_merged[nx - k0, 0]) fft_merged[0, k0] = np.conj(fft_merged[0, nx - k0]) - freq_i = np.tile(np.fft.fftfreq(array_high.shape[0], d=1 / ds_factor), nx).reshape( - (nx, nx) - ) + # Compute frequency arrays + freq_i = np.fft.fftfreq(nx, d=1 / ds_factor) + freq_i = np.tile(freq_i, (nx, 1)) freq_j = freq_i.T + # Compute frequency domain adjustment ddx = np.pi * (1 / nax - 1 / nx) / np.abs(freq_i[0, 1] - freq_i[0, 0]) - fx2 = freq_array_high**2 - fax2 = freq_array_low[k0, k0] ** 2 + freq_squared_high = freq_array_high**2 + freq_squared_low_center = freq_array_low[k0, k0] ** 2 - fft_merged = fft_array_high * (fx2 > fax2) + fft_merged * (fx2 <= fax2) * np.exp( + # Fuse in the frequency domain + mask_high = freq_squared_high > freq_squared_low_center + mask_low = ~mask_high + fft_merged = fft_high * mask_high + fft_merged * mask_low * np.exp( -1j * ddx * freq_i - 1j * ddx * freq_j ) + # Inverse FFT to obtain the merged array in the spatial domain merged = np.real(np.fft.ifftn(fft_merged)) / len(fft_merged) return merged @@ -143,6 +190,10 @@ def _make_gaussian_kernel(ds_factor): def _balanced_spatial_average(array, kernel): + """ + Compute the balanced spatial average of an array using a given kernel while handling + missing or invalid values. + """ array = array.copy() mask_valid = np.isfinite(array) array[~mask_valid] = 0.0 @@ -214,62 +265,64 @@ def downscale( """ - if spectral_fusion: - precip_transformed = _gaussianize(precip) - else: - precip_transformed = precip + # Validate inputs + if not np.isfinite(precip).all(): + raise ValueError("All values in 'precip' must be finite.") + if not isinstance(ds_factor, int) or ds_factor <= 0: + raise ValueError("'ds_factor' must be a positive integer.") + + # Preprocess the input field if spectral fusion is enabled + precip_transformed = _gaussianize(precip) if spectral_fusion else precip + # Compute frequency arrays for the original and high-resolution fields freq_array = _compute_freq_array(precip_transformed) freq_array_highres = _compute_freq_array(precip_transformed, ds_factor) + # Estimate spectral slope alpha if not provided if alpha is None: alpha = _estimate_alpha(precip_transformed, freq_array) - white_noise_field = np.random.rand(*freq_array_highres.shape) - white_noise_field_complex = np.exp(complex(0, 1) * 2 * np.pi * white_noise_field) - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - noise_field_complex = white_noise_field_complex * np.sqrt( - freq_array_highres**-alpha - ) - noise_field_complex[0, 0] = 0 - noise_field = np.fft.ifft2(noise_field_complex).real + # Generate noise field + noise_field = _compute_noise_field(freq_array_highres, alpha) + # Apply spectral fusion if enabled if spectral_fusion: noise_field /= noise_field.shape[0] ** 2 noise_field = np.exp(noise_field) - noise_field = _apply_spectral_fusion( precip_transformed, noise_field, freq_array, freq_array_highres, ds_factor ) + # Normalize and exponentiate the noise field noise_field /= noise_field.std() noise_field = np.exp(noise_field) + # Aggregate the noise field to low resolution noise_lowres = aggregate_fields(noise_field, ds_factor, axis=(0, 1)) - precip = np.kron(precip, np.ones((ds_factor, ds_factor))) - noise_lowres = np.kron(noise_lowres, np.ones((ds_factor, ds_factor))) + # Expand input and noise fields to high resolution + precip_expanded = np.kron(precip, np.ones((ds_factor, ds_factor))) + noise_lowres_expanded = np.kron(noise_lowres, np.ones((ds_factor, ds_factor))) - if kernel_type is not None: - try: - kernel = _make_kernel[kernel_type](ds_factor) - except KeyError: + # Apply smoothing if a kernel type is provided + if kernel_type: + if kernel_type not in _make_kernel: raise ValueError( - f"kernel type '{kernel_type}' is invalid, " - f"available kernels: {list(_make_kernel)}" + f"kernel type '{kernel_type}' is invalid, available kernels: {list(_make_kernel)}" ) + kernel = _make_kernel[kernel_type](ds_factor) + precip_expanded = _balanced_spatial_average(precip_expanded, kernel) + noise_lowres_expanded = _balanced_spatial_average(noise_lowres_expanded, kernel) - precip = _balanced_spatial_average(precip, kernel) - noise_lowres = _balanced_spatial_average(noise_lowres, kernel) - - norm_k0 = precip / noise_lowres - + # Normalize the high-res precipitation field by the low-res noise field + norm_k0 = precip_expanded / noise_lowres_expanded precip_highres = noise_field * norm_k0 + # Apply thresholding if specified if threshold is not None: precip_highres[precip_highres < threshold] = 0 + # Return the downscaled field and optionally the spectral slope alpha if return_alpha: return precip_highres, alpha From 88275153371a6354009daf539c57ede7e13982a8 Mon Sep 17 00:00:00 2001 From: ecasellas Date: Mon, 10 Jun 2024 13:23:35 +0200 Subject: [PATCH 21/21] Correction of merged normalisation --- pysteps/downscaling/rainfarm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py index d9906390d..b40e359df 100644 --- a/pysteps/downscaling/rainfarm.py +++ b/pysteps/downscaling/rainfarm.py @@ -157,7 +157,7 @@ def compute_psd(array, fft_size): ) # Inverse FFT to obtain the merged array in the spatial domain - merged = np.real(np.fft.ifftn(fft_merged)) / len(fft_merged) + merged = np.real(np.fft.ifftn(fft_merged)) / fft_merged.size return merged