fixed inconsistent return types; added testing

Author: btschwertfeger

README.md

@@ -1,31 +1,43 @@
 # Bias-Adjustment-Python
 
-[![Generic badge](https://img.shields.io/badge/license-MIT-green.svg)](https://shields.io/)
-[![Generic badge](https://img.shields.io/badge/python-3.7+-blue.svg)](https://shields.io/)
+<div style="text-align: center">
+
 [![GitHub](https://badgen.net/badge/icon/github?icon=github&label)](https://github.com/btschwertfeger/Bias-Adjustment-Python)
+[![Generic badge](https://img.shields.io/badge/python-3.7+-blue.svg)](https://shields.io/)
+[![Generic badge](https://img.shields.io/badge/license-MIT-green.svg)](https://shields.io/)
+[![Downloads](https://pepy.tech/badge/python-cmethods)](https://pepy.tech/project/python-cmethods)
+[![Downloads](https://pepy.tech/badge/python-cmethods/month)](https://pepy.tech/project/python-cmethods)
 
+</div>
 Collection of different scale- and distribution-based bias adjustment techniques for climatic research. (see `examples.ipynb` for help)
 
 Bias adjustment procedures in Python are very slow, so they should not be used on large data sets.
 A C++ implementation that works way faster can be found here: [https://github.com/btschwertfeger/Bias-Adjustment-Cpp](https://github.com/btschwertfeger/Bias-Adjustment-Cpp).
-____
+
+---
+
 ## Available methods:
+
 - Linear Scaling (additive and multiplicative)
 - Variance Scaling (additive)
 - Delta (Change) Method (additive and multiplicative)
 - Quantile Mapping (additive)
 - Detrended Quantile Mapping (additive and multiplicative)
 - Quantile Delta Mapping (additive and multuplicative)
 
-____
+---
+
 ## Usage
 
 ### Installation
+
 ```bash
 python3 -m pip install python-cmethods
 ```
+
 ### Import and application
-```python 
+
+```python
 import xarray as xr
 from cmethods.CMethods import CMethods
 cm = CMethods()
@@ -41,25 +53,27 @@ ls_result = cm.linear_scaling(
     kind = '+' # *
 )
 
-qdm_result = cm.adjust_2d( # 2d = 2 spatial and 1 time dimension
+qdm_result = cm.adjust_3d( # 3d = 2 spatial and 1 time dimension
     method = 'quantile_delta_mapping',
     obs = obsh['tas'],
     simh = simh['tas'],
     simp = simp['tas'],
     n_quaniles = 1000,
     kind = '+' # *
 )
-# * to calculate the relative rather than the absolute change, 
+# * to calculate the relative rather than the absolute change,
 # '*' can be used instead of '+' (this is prefered when adjusting
 # ratio based variables like precipitation)
 ```
 
-____
+---
+
 ## Examples (see repository on [GitHub](https://github.com/btschwertfeger/Bias-Adjustment-Python))
 
 `/examples/examples.ipynb`: Notebook containing different methods and plots
 
 `/examples/do_bias_correction.py`: Example script for adjusting climate data
+
 ```bash
 python3 do_bias_correction.py   \
     --obs input_data/obs.nc     \
@@ -69,17 +83,21 @@ python3 do_bias_correction.py   \
     --variable tas              \
     --unit '°C'                 \
     --group time.month          \
-    --kind + 
+    --kind +
 ```
 
 - Linear and variance, as well as delta change method require `--group time.month` as argument.
 - Adjustment methods that apply changes in distributional biasses (QM, QDM, DQM; EQM, ...) need the `--nquantiles` argument set to some integer.
 - Data sets should have the same spatial resolutions.
-____
+
+---
+
 ## Notes:
+
 - Computation in Python takes some time, so this is only for demonstration. When adjusting large datasets, its best to the C++ implementation mentioned above.
 
 ## References
+
 - Schwertfeger, Benjamin Thomas (2022) The influence of bias corrections on variability, distribution, and correlation of temperatures in comparison to observed and modeled climate data in Europe (https://epic.awi.de/id/eprint/56689/)
 - Linear Scaling and Variance Scaling based on: Teutschbein, Claudia and Seibert, Jan (2012) Bias correction of regional climate model simulations for hydrological climate-change impact studies: Review and evaluation of different methods (https://doi.org/10.1016/j.jhydrol.2012.05.052)
 - Delta Method based on: Beyer, R. and Krapp, M. and Manica, A.: An empirical evaluation of bias correction methods for palaeoclimate simulations (https://doi.org/10.5194/cp-16-1493-2020)

cmethods/CMethods.py

@@ -23,7 +23,7 @@
     simp = data to correct (predicted simulated data) ($T_{sim,p}$)                     \
                                                                                         \
                                                                                         \
-    F = Cummulative Distribution Function                                               \
+    F = Cumulative Distribution Function                                               \
     \mu = mean                                                                          \
     \sigma = standard deviation                                                         \
     i = index                                                                           \
@@ -46,9 +46,14 @@ def __init__(self, method: str, available_methods: list):
 
     SCALING_METHODS = ['linear_scaling', 'variance_scaling', 'delta_method']
     DISTRIBUTION_METHODS = ['quantile_mapping', 'quantile_delta_mapping']
+    
     CUSTOM_METHODS = SCALING_METHODS + DISTRIBUTION_METHODS
+    
     METHODS = CUSTOM_METHODS #+ XCLIM_SDBA_METHODS
 
+    ADDITIVE = ['+', 'add']
+    MULTIPLICATIVE = ['*', 'mult']
+
     def __init__(self):
         pass
 
@@ -68,7 +73,7 @@ def get_function(cls, method: str):
         else: raise UnknownMethodError(method, cls.METHODS)
 
     @classmethod
-    def adjust_2d(cls,
+    def adjust_3d(cls,
         method: str,
         obs: xr.core.dataarray.DataArray,
         simh: xr.core.dataarray.DataArray,
@@ -104,9 +109,9 @@ def adjust_2d(cls,
             > simh = xarray.open_dataset('path/to/simulated/data.nc')
             > obs = xarray.open_dataset('path/to/observed/data.nc')
             > simp = xarray.open_dataset('path/to/simulated_future/data.nc')
-            > variable = 'temperature'
+            > variable = 'tas'
 
-            > adjusted_data = CMethods().adjust_2d(
+            > adjusted_data = CMethods().adjust_3d(
                 method = 'quantile_delta_mapping',
                 obs = obs[variable],
                 simh = simh[variable],
@@ -121,7 +126,7 @@ def adjust_2d(cls,
         simh = simh.transpose('lat', 'lon', 'time')
         simp = simp.transpose('lat', 'lon', 'time')
 
-        if group == None and method in SCALING_METHODS: group = 'time.month'
+        if group == None and method in cls.SCALING_METHODS: group = 'time.month'
 
         result = simp.copy(deep=True).load()
         len_lat, len_lon = len(obs.lat), len(obs.lon)
@@ -160,7 +165,7 @@ def adjust_2d(cls,
     @classmethod
     def pool_adjust(cls, params) -> xr.core.dataarray.DataArray:
         ''' Adjustment along longitude for one specific latitude
-            used by cls.adjust_2d as callbackfunction for multiprocessing.Pool
+            used by cls.adjust_3d as callbackfunction for multiprocessing.Pool
         '''
 
         method = params['method']
@@ -249,7 +254,7 @@ def linear_scaling(cls,
             simh (xarray.core.dataarray.DataArray): simulated historical Data
             simp (xarray.core.dataarray.DataArray): future simulated Data
             group (str): [optional] Group / Period (e.g.: 'time.month')
-            kind (str): '+' or '*', default: '+'
+            kind (str): [optional] '+' or '*', default: '+'
 
         ----- R E T U R N -----
 
@@ -259,13 +264,13 @@ def linear_scaling(cls,
             > obs = xarray.open_dataset('path/to/observed/data.nc')
             > simh = xarray.open_dataset('path/to/simulated/data.nc')
             > simp = xarray.open_dataset('path/to/predicted/data.nc')
-            > variable = 'temperature'
+            > variable = 'tas'
 
             > result = CMethods().linear_scaling(
             >    obs=obs[variable],
             >    simh=simh[variable],
             >    simp=simp[variable],
-            >    group='time.month'
+            >    group='time.month' # optional
             >)
 
         ----- E Q U A T I O N S -----
@@ -281,11 +286,10 @@ def linear_scaling(cls,
             https://doi.org/10.1016/j.jhydrol.2012.05.052
 
         '''
-
         if group != None: return cls.grouped_correction(method='linear_scaling', obs=obs, simh=simh, simp=simp, group=group, kind=kind, **kwargs)
         else:
-            if kind == '+': return np.array(simp) + (np.nanmean(obs) - np.nanmean(simh)) # Eq. 1
-            elif kind == '*': return np.array(simp) * (np.nanmean(obs) / np.nanmean(simh)) # Eq. 2
+            if kind in cls.ADDITIVE: return np.array(simp) + (np.nanmean(obs) - np.nanmean(simh)) # Eq. 1
+            elif kind in cls.MULTIPLICATIVE: return np.array(simp) * (np.nanmean(obs) / np.nanmean(simh)) # Eq. 2
             else: raise ValueError('Scaling type invalid. Valid options for param kind: "+" and "*"')
 
     # ? -----========= V A R I A N C E - S C A L I N G =========------
@@ -306,7 +310,7 @@ def variance_scaling(cls,
             simh (xarray.core.dataarray.DataArray): simulated historical Data
             simp (xarray.core.dataarray.DataArray): future simulated Data
             group (str): [optional] Group / Period (e.g.: 'time.month')
-            kind (str): '+' or '*', default: '+'
+            kind (str): '+' or '*', default: '+' # '*' is not implemented
 
         ----- R E T U R N -----
 
@@ -316,7 +320,7 @@ def variance_scaling(cls,
             > obs = xarray.open_dataset('path/to/observed/data.nc')
             > simh = xarray.open_dataset('path/to/simulated/data.nc')
             > simp = xarray.open_dataset('path/to/predicted/data.nc')
-            > variable = 'temperature'
+            > variable = 'tas'
 
             > result = CMethods().variance_scaling(obs=obs[variable], simh=simh[variable], simp=simp[variable] group='time.dayofyear')
 
@@ -341,8 +345,8 @@ def variance_scaling(cls,
         '''
         if group != None: return cls.grouped_correction(method='variance_scaling', obs=obs, simh=simh, simp=simp, group=group, kind=kind, **kwargs)
         else:
-            LS_simh = cls.linear_scaling(obs, simh, simh)               # Eq. 1
-            LS_simp = cls.linear_scaling(obs, simh, simp)               # Eq. 2
+            LS_simh = cls.linear_scaling(obs, simh, simh, group=None)   # Eq. 1
+            LS_simp = cls.linear_scaling(obs, simh, simp, group=None)   # Eq. 2
 
             VS_1_simh = LS_simh - np.nanmean(LS_simh)                   # Eq. 3
             VS_1_simp = LS_simp - np.nanmean(LS_simp)                   # Eq. 4
@@ -379,7 +383,7 @@ def delta_method(cls,
             > obs = xarray.open_dataset('path/to/observed/data.nc')
             > simh = xarray.open_dataset('path/to/simulated/data.nc')
             > simp = xarray.open_dataset('path/to/predicted/data.nc')
-            > variable = 'temperature'
+            > variable = 'tas'
 
             > result = CMethods().delta_method(obs=obs[variable], simh=simh[variable], group='time.month')
 
@@ -399,8 +403,8 @@ def delta_method(cls,
         '''
         if group != None: return cls.grouped_correction(method='delta_method', obs=obs, simh=simh, simp=simp, group=group, kind=kind, **kwargs)
         else:
-            if kind == '+': return np.array(obs) + (np.nanmean(simp) - np.nanmean(simh))     # Eq. 1
-            elif kind == '*': return np.array(obs) * (np.nanmean(simp) / np.nanmean(simh))   # Eq. 2
+            if kind in cls.ADDITIVE: return np.array(obs) + (np.nanmean(simp) - np.nanmean(simh))     # Eq. 1
+            elif kind in cls.MULTIPLICATIVE: return np.array(obs) * (np.nanmean(simp) / np.nanmean(simh))   # Eq. 2
             else: raise ValueError(f'{kind} not implemented! Use "+" or "*" instead.')
 
 
@@ -468,7 +472,7 @@ def quantile_mapping(cls,
         cdf_obs = cls.get_cdf(obs, xbins)
         cdf_simh = cls.get_cdf(simh, xbins)
 
-        if kwargs.get('detrended', False) or kind in [ '*', 'mult' ]:
+        if kwargs.get('detrended', False) or kind in cls.MULTIPLICATIVE:
             '''detrended => shift mean of $X_{sim,p}$ to range of $X_{sim,h}$ to adjust extremes'''
             for month, idxs in res.groupby('time.month').groups.items():
                 m_simh, m_simp = [], []
@@ -481,7 +485,7 @@ def quantile_mapping(cls,
                 m_simh_mean = np.nanmean(m_simh)     
                 m_simp_mean = np.nanmean(m_simp) 
 
-                if kind == '+':
+                if kind in cls.ADDITIVE:
                     epsilon = np.interp(m_simp - m_simp_mean + m_simh_mean, xbins, cdf_simh)         # Eq. 1     
                     X = cls.get_inverse_of_cdf(cdf_obs, epsilon, xbins) + m_simp_mean - m_simh_mean   # Eq. 1     
                 else: 
@@ -490,7 +494,7 @@ def quantile_mapping(cls,
                     #x = cm.get_inverse_of_cdf(cdf_obs, epsilon, xbins) * (m_simp_mean / m_simh_mean)
                 for i, idx in enumerate(idxs): res.values[idx] = X[i]
             return res
-        elif kind in [ '+', 'add' ]: # additive, no detrend 
+        elif kind in cls.ADDITIVE: # additive, no detrend 
             epsilon = np.interp(simp, xbins, cdf_simh)                                 # Eq. 1            
             res.values = cls.get_inverse_of_cdf(cdf_obs, epsilon, xbins)               # Eq. 1      
             return res
@@ -527,7 +531,7 @@ def empirical_quantile_mapping(cls,
             > obs = xarray.open_dataset('path/to/observed/data.nc')
             > simh = xarray.open_dataset('path/to/simulated/data.nc')
             > simp = xarray.open_dataset('path/to/future/data.nc')
-            > variable = 'temperature'
+            > variable = 'tas'
             > result = CMethods().empirical_quantile_mapping(
             >    obs=obs[variable],
             >    simh=simh[variable],
@@ -613,7 +617,8 @@ def quantile_delta_mapping(cls,
         '''
 
         if group != None: return cls.grouped_correction(method='quantile_delta_mapping', obs=obs, simh=simh, simp=simp, group=group, n_quantiles=n_quantiles, kind=kind, **kwargs)
-        elif kind == '+':
+        elif kind in cls.ADDITIVE:
+            res = simp.copy(deep=True)
             obs, simh, simp = np.array(obs), np.array(simh), np.array(simp) # to achieve higher accuracy
             global_max = max(np.amax(obs), np.amax(simh))
             global_min = min(np.amin(obs), np.amin(simh))
@@ -628,9 +633,11 @@ def quantile_delta_mapping(cls,
             epsilon = np.interp(simp, xbins, cdf_simp)                            # Eq. 1.1
             QDM1 = cls.get_inverse_of_cdf(cdf_obs, epsilon, xbins)                # Eq. 1.2
             delta = simp - cls.get_inverse_of_cdf(cdf_simh, epsilon, xbins)       # Eq. 1.3
-            return QDM1 + delta                                                   # Eq. 1.4
+            res.values = QDM1 + delta                                             # Eq. 1.4
+            return res
 
-        elif kind == '*':
+        elif kind in cls.MULTIPLICATIVE:
+            res = simp.copy(deep=True)
             obs, simh, simp = np.array(obs), np.array(simh), np.array(simp)
             global_max = max(np.amax(obs), np.amax(simh))
             wide = global_max / n_quantiles
@@ -643,8 +650,8 @@ def quantile_delta_mapping(cls,
             epsilon = np.interp(simp, xbins, cdf_simp)                            # Eq. 1.1
             QDM1 = cls.get_inverse_of_cdf(cdf_obs, epsilon, xbins)                # Eq. 1.2
             delta = simp / cls.get_inverse_of_cdf(cdf_simh, epsilon, xbins)       # Eq. 2.3
-            return QDM1 * delta                                                   # Eq. 2.4
-
+            res.values = QDM1 * delta                                             # Eq. 2.4
+            return res
         else: raise ValueError(f'Unknown kind {kind}!')
 
     # * -----========= G E N E R A L  =========------

cmethods/__version__.py

@@ -1,2 +1,2 @@
-VERSION = (0,5,1)
+VERSION = (0,5,3)
 __version__ = '.'.join(map(str, VERSION))