removed xclim functions; added setup.py; minor change in project structure

Author: btschwertfeger

.gitignore

@@ -1,2 +1,7 @@
 .DS_Store
-__pycache__
+__pycache__
+.ipynb_checkpoints
+
+build/
+dist/
+*.egg-info/

MANIFEST.in

@@ -0,0 +1 @@
+include README.md LICENSE

README.md

@@ -1,10 +1,9 @@
 # Bias-Adjustment-Python
 
-Collection of different scale- and distribution-based bias adjustments for climatic research. This methods are part of the bachelor thesis of Benjamin T. Schwertfeger.
-During this thesis, many of these methods have also been implemented in C++.
-This can be found here: [https://github.com/btschwertfeger/Bias-Adjustment-Cpp](https://github.com/btschwertfeger/Bias-Adjustment-Cpp).
+Collection of different scale- and distribution-based bias adjustment techniques for climatic research. (see `examples.ipynb` for help)
 
-There is also a Jupyter Notebook that serves as example.
+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).
 ____
 ## Run adjustment:
 ```bash
@@ -28,14 +27,6 @@ ____
 |Quantile Mapping|quantile_mapping|
 |Quantile Delta Mapping|quantile_delta_mapping|
 
-## Methods adapted from [xclim](https://xclim.readthedocs.io/en/stable/sdba.html):
-|Method| `--method` parameter|
-|-----|-----|
-|Empirical Quantile Mapping|xclim_eqm|
-|Detrended Quantile Mapping|xclim_dqm|
-|Quantile Delta Mapping|xclime_qdm|
-
-
 ____
 # Notes:
 - Linear and variance, as well as delta change method require `--group time.month` as argument.

CMethods.py renamed to cmethods/CMethods.py

@@ -1,17 +1,11 @@
 #!/bin/python3
 
 import multiprocessing
-import xclim as xc
-from xclim import sdba
 import xarray as xr
 import numpy as np
-from scipy.stats.mstats import mquantiles
-from scipy.interpolate import interp1d
-
 from tqdm import tqdm
 
-
-__descrption__ = 'Script to adjust climate biases in 3D Climate data'
+__descrption__ = 'Script to adjust climate bias in climate data'
 __author__ = 'Benjamin Thomas Schwertfeger'
 __copyright__ = __author__
 __email__ = '[email protected]'
@@ -52,10 +46,9 @@ def __init__(self, method: str, available_methods: list):
 
     CUSTOM_METHODS = [
         'linear_scaling', 'variance_scaling', 'delta_method', 
-        'quantile_mapping', 'empirical_quantile_mapping', 'quantile_delta_mapping'
+        'quantile_mapping', 'quantile_delta_mapping'
     ]
-    XCLIM_SDBA_METHODS = ['xclim_eqm', 'xclim_dqm', 'xclim_qdm']
-    METHODS = CUSTOM_METHODS + XCLIM_SDBA_METHODS
+    METHODS = CUSTOM_METHODS #+ XCLIM_SDBA_METHODS
 
     def __init__(self):
         pass
@@ -67,10 +60,7 @@ def get_available_methods(cls) -> list:
 
     @classmethod
     def get_function(cls, method: str):
-        if method == 'xclim_dqm': return sdba.adjustment.DetrendedQuantileMapping.train 
-        elif method == 'xclim_eqm': return sdba.adjustment.EmpiricalQuantileMapping.train        
-        elif method == 'xclim_qdm': return sdba.adjustment.QuantileDeltaMapping.train 
-        elif method == 'linear_scaling': return cls.linear_scaling        
+        if method == 'linear_scaling': return cls.linear_scaling        
         elif method == 'variance_scaling': return cls.variance_scaling        
         elif method == 'delta_method':  return cls.delta_method
         elif method == 'quantile_mapping': return cls.quantile_mapping
@@ -140,40 +130,7 @@ def adjust_2d(cls,
         result = simp.copy(deep=True).load()
         len_lat, len_lon = len(obs.lat), len(obs.lon)
 
-        if method in cls.XCLIM_SDBA_METHODS:
-            if n_jobs == 1:
-                method = cls.get_function(method)   
-                for lat in tqdm(range(len_lat)):
-                    for lon in range(len_lon):
-                        result[lat,lon] = method(
-                            obs = obs[lat,lon], 
-                            simh = simh[lat,lon], 
-                            simp = simp[lat,lon],
-                            tslice_adjust = tslice_adjust, 
-                            n_quantiles = n_quantiles, 
-                            kind = kind, 
-                            group = group, 
-                            window = window,
-                            **kwargs
-                        ) 
-            else:
-                with multiprocessing.Pool(processes=n_jobs) as pool:
-                    params: [dict] = [{
-                        'method': method, 
-                        'obs': obs[lat], 
-                        'simh': simh[lat], 
-                        'simp': simp[lat], 
-                        'tslice_adjust': tslice_adjust, 
-                        'kind': kind, 
-                        'group': group, 
-                        'window': window,
-                        'kwargs': kwargs
-                    } for lat in range(len_lat)]
-                    for lat, corrected in enumerate(pool.map(cls.pool_adjust, params)):
-                        result[lat] = corrected
-            return result.transpose('time', 'lat', 'lon')
-        
-        elif method in cls.CUSTOM_METHODS:
+        if method in cls.CUSTOM_METHODS:
             if n_jobs == 1:
                 method = cls.get_function(method)
                 for lat in tqdm(range(len_lat)):
@@ -226,25 +183,7 @@ def pool_adjust(cls, params) -> xr.core.dataarray.DataArray:
         len_lon = len(obs.lon)
 
         func_adjustment = None
-        if method in cls.XCLIM_SDBA_METHODS:
-            func_adjustment = cls.get_function(method)
-            for lon in range(len_lon):
-                result[lon] = cls.xclim_sdba_adjustment(
-                    method = func_adjustment, 
-                    method_name = method, 
-                    obs = obs, 
-                    simh = simh, 
-                    simp = simp, 
-                    tslice_adjust = tslice_adjust,
-                    n_quantiles = n_quantiles, 
-                    group = group, 
-                    window = window,
-                    save_model = save_model,
-                    **kwargs
-                )
-            return result
-            
-        elif method in cls.CUSTOM_METHODS:
+        if method in cls.CUSTOM_METHODS:
             func_adjustment = cls.get_function(method)                
             kwargs['n_quantiles'] = n_quantiles
             kwargs['kind'] = kind
@@ -255,57 +194,7 @@ def pool_adjust(cls, params) -> xr.core.dataarray.DataArray:
             return result
 
         else: raise UnknownMethodError(method, cls.METHODS)       
-            
-    @staticmethod
-    def xclim_sdba_adjustment(
-        method,
-        method_name: str,
-        obs: xr.core.dataarray.DataArray, 
-        simh: xr.core.dataarray.DataArray, 
-        simp: xr.core.dataarray.DataArray,
-        tslice_adjust: slice=None,
-        n_quantiles: int=100, 
-        kind: str='+', 
-        group: str='time.month', 
-        window: int=1,
-        save_model: bool=False
-    ) -> xr.core.dataarray.DataArray:
-        '''Method to adjust 1 dimensional climate data using the xclim.sdba library
-        
-        Note: This Method should only be called by the method adjust_2d.
-        
-        ----- P A R A M E T E R S -----
-        
-            method (method): adjustment method (once of the xclim.sdba library)
-            method_name (str): Name of the method to use in filename of the model
-    
-            obs (xarray.core.dataarray.DataArray): observed / obserence Data
-            simh (xarray.core.dataarray.DataArray): simulated historical Data
-            simp (xarray.core.dataarray.DataArray): future simulated Data to adjust
-            
-            tslice_adjust (slice): Timespan to adjust, default: None
-            n_quantiles (int): Number of quantiles to involve
-            kind (str): Kind of adjustment ('+' or '*'), default: '+'
-            group (str): Group data by, default: 'time.month'
-            window (int): Grouping window, default: 1
-            
-        ----- R E T U R N -----
-            
-            xarray.core.dataarray.DataArray: Adjusted data 
-                  
-        '''
-        
-        grouper = xc.sdba.Grouper(group=group, window=window) 
-        model = method(obs, simh, nquantiles=n_quantiles, kind=kind, group=grouper)
-        
-        if save_model: model.ds.to_netcdf(f'{method_name}_model_nquantiles-{n_quantiles}_kind-{kind}_group-{group}_window-{window}.nc')
-        if tslice_adjust != None: simp = simp.sel(time=tslice_adjust)
-
-        bias = model.adjust(simp)
-        bias.attrs = simp.attrs
-        
-        return bias 
-    
+                
     @classmethod 
     def grouped_correction(cls,
         method: str,
@@ -458,7 +347,6 @@ 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:
-            obs, simh, simp = np.array(obs), np.array(simh), np.array(simp)
             LS_simh = cls.linear_scaling(obs, simh, simh)               # Eq. 1
             LS_simp = cls.linear_scaling(obs, simh, simp)               # Eq. 2
 
@@ -544,19 +432,20 @@ def quantile_mapping(cls,
 
         if group != None: return cls.grouped_correction(method='quantile_mapping', obs=obs, simh=simh, simp=simp, group=group, n_quantiles=n_quantiles, kind=kind, **kwargs)
         elif kind == '+':
-            # make np.array to achieve higher accuracy (idk why)
+            # create np.array to achieve higher accuracy (idk why)
             obs, simh, simp = np.array(obs), np.array(simh), np.array(simp)
-            # define quantile boundaries (xbins)
+            
             global_max = max(np.amax(obs), np.amax(simh))
-            global_min = min(np.amin(obs), np.amin(simh)) # change to 0.0 if precipitation
-            wide = abs(global_max - global_min) / n_quantiles # change to global_max/n_quantiles if precipitation
+            global_min = min(np.amin(obs), np.amin(simh)) 
+            wide = abs(global_max - global_min) / n_quantiles 
             xbins = np.arange(global_min, global_max + wide, wide)
 
             cdf_obs = cls.get_cdf(obs, xbins)
             cdf_simh = cls.get_cdf(simh, xbins)
             epsilon = np.interp(simp, xbins, cdf_simh)                                 # Eq. 1
 
             return cls.get_inverse_of_cdf(cdf_obs, epsilon, xbins)                     # Eq. 1
+            
         elif kind == '*': 
             ''' Inspired by Adrian Tompkins [email protected] posted here: 
                 https://www.researchgate.net/post/Does-anyone-know-about-bias-correction-and-quantile-mapping-in-PYTHON
@@ -625,7 +514,7 @@ def empirical_quantile_mapping(cls,
             https://svn.oss.deltares.nl/repos/openearthtools/trunk/python/applications/hydrotools/hydrotools/statistics/bias_correction.py
             
         '''
-        raise ValueError('idk if it is allowed to use this so please have a look at:https://svn.oss.deltares.nl/repos/openearthtools/trunk/python/applications/hydrotools/hydrotools/statistics/bias_correction.py ')
+        raise ValueError('idk if it is allowed to use this so please have a look at: https://svn.oss.deltares.nl/repos/openearthtools/trunk/python/applications/hydrotools/hydrotools/statistics/bias_correction.py ')
         # if group != None:
         #     return cls.grouped_correction(
         #         method = 'empirical_quantile_mapping', 
@@ -689,7 +578,6 @@ def quantile_delta_mapping(cls,
             # calculate exact cdf values of $F_{sim,p}[T_{sim,p}(t)]$
             epsilon = np.interp(simp, xbins, cdf_simp)                            # Eq. 1
 
-            # invert F_{obs,h}
             QDM1 = cls.get_inverse_of_cdf(cdf_obs, epsilon, xbins)                # Eq. 2
             delta = simp - cls.get_inverse_of_cdf(cdf_simh, epsilon, xbins)       # Eq. 3 
 
@@ -712,7 +600,7 @@ def get_cdf(a, xbins: list) -> np.array:
 
     @staticmethod    
     def get_inverse_of_cdf(base_cdf, insert_cdf, xbins) -> np.array:
-        ''' returns the inverse cummulative distribution function of base_cdf ($$F_{base_cdf}\left[insert_cdf\right])$$'''
+        ''' returns the inverse cummulative distribution function of base_cdf ($F_{base_cdf}\left[insert_cdf\right])$'''
         return np.interp(insert_cdf, base_cdf, xbins)
 
     @staticmethod

cmethods/__init__.py

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

cmethods/__version__.py

@@ -0,0 +1,4 @@
+xarray
+numpy
+tqdm
+matplotlib

cmethods/requirements.py

@@ -11,7 +11,7 @@
 import logging, sys
 import xarray as xr
 
-from CMethods import CMethods
+from cmethods.CMethods import CMethods
 
 # * ----- L O G G I N G -----
 formatter = logging.Formatter(