sadpy

Python implementation of the "scale, anisotropy, direction" approach to spatial forecast verification based on two-dimensional wavelet transforms (https://doi.org/10.1002/qj.3964).

Installation

git clone https://github.com/s6sebusc/sadpy.git
cd sadpy
pip install .

Usage example

Suppose you have 2D arrays forecast and observation, then you could say

from sadpy import wv_analysis, sad

obs_analysis = wv_analysis(observation)
for_analysis = wv_analysis(forecast)

scores = sad( for_analysis, obs_analysis )

The result in scores is a dict of the scale, anisotropy and direction errors.

Provenance

The original code for this method is available in the R-package sad (https://cran.r-project.org/package=sad). To avoid dependencies, this package also contains a limited implementation of the dual-tree complex wavelet transform (dtcwt, https://doi.org/10.1109/MSP.2005.1550194). It is not copied from the dtcwt library (https://github.com/xir4n/dtcwt) but based on the dualtrees R-package (https://CRAN.R-project.org/package=dualtrees) which originally took inspiration from dtcwt.

Alternative wavelet transforms

The sadpy.dtcwt implementation can be replaced by any of the more mature and powerful dual-tree implementations out there (https://github.com/xir4n/dtcwt, https://github.com/fbcotter/pytorch_wavelets, ...). It can also be replaced by any other wavelet transform, for example from pywavelets. Anisotropy and direction are, however, only available with the dtcwt.

To replace sadpy.dtcwt by a discrete haar wavelet transform from pywavelets, you could say

from pywt import wavedec2

def dwt(x, n_levels, wavelet='haar'):
    res = wavedec2(x, wavelet, level=n_levels, axes=(-2,-1))
    res = [ np.stack( res[n_levels-i], axis=-1 ) for i in range(n_levels) ]
    return res

obs_analysis = wv_analysis(observation, wv_transform=dwt)