fixed detrended quantile mapping for ratio based variables

Author: btschwertfeger

README.md

@@ -8,13 +8,15 @@ Collection of different scale- and distribution-based bias adjustment techniques
 
 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)
 
-## Available methods
-- Linear Scaling
-- Variance Scaling 
-- Delta (Change) Method
-- Quantile and Detrended Quantile Mapping
-- Quantile Delta Mapping
 ____
 ## Usage
 
@@ -24,6 +26,7 @@ python3 -m pip install python-cmethods
 ```
 ### Import and application
 ```python 
+import xarray as xr
 from cmethods.CMethods import CMethods
 cm = CMethods()
 

cmethods/CMethods.py

@@ -44,10 +44,9 @@ class UnknownMethodError(Exception):
         def __init__(self, method: str, available_methods: list):
             super().__init__(f'Unknown method "{method}"! Available methods: {available_methods}')
 
-    CUSTOM_METHODS = [
-        'linear_scaling', 'variance_scaling', 'delta_method',
-        'quantile_mapping', 'quantile_delta_mapping'
-    ]
+    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
 
     def __init__(self):
@@ -122,6 +121,8 @@ 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'
+
         result = simp.copy(deep=True).load()
         len_lat, len_lon = len(obs.lat), len(obs.lon)
 
@@ -169,7 +170,6 @@ def pool_adjust(cls, params) -> xr.core.dataarray.DataArray:
         n_quantiles = params.get('n_quantiles', 100)
         kind = params.get('kind', '+')
         group = params.get('group', None)
-        save_model = params.get('save_model', None)
         kwargs = params.get('kwargs', {})
 
         result = simp.copy(deep=True).load()
@@ -393,7 +393,7 @@ def delta_method(cls,
         ----- R E F E R E N C E S -----
             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)
 
-            and 
+            and
             https://svn.oss.deltares.nl/repos/openearthtools/trunk/python/applications/hydrotools/hydrotools/statistics/bias_correction.py
 
         '''
@@ -435,67 +435,65 @@ def quantile_mapping(cls,
             Add (+):
                 (1.) X^{*QM}_{sim,p}(i) = F^{-1}_{obs,h} \left\{F_{sim,h}\left[X_{sim,p}(i)\right]\right\}
             Mult (*):
-                (1.) --//--
+                (2.) X^{*QM}_{sim,p}(i) = F^{-1}_{obs,h}
+                        \Biggl\{
+                            F_{sim,h}\left[
+                                \frac{
+                                    \mu{X_{sim,h}} \mu{X_{sim,p}(i)}
+                                }{
+                                    \mu{X_{sim,p}(i)}
+                                }
+                            \right]
+                        \Biggr\}
+                        \frac{
+                            \mu{X_{sim,p}(i)}
+                        }{
+                            \mu{X_{sim,h}}
+                        }  
 
         ----- R E F E R E N C E S -----
             Alex J. Cannon and Stephen R. Sobie and Trevor Q. Murdock Bias Correction of GCM Precipitation by Quantile Mapping: How Well Do Methods Preserve Changes in Quantiles and Extremes?
             https://doi.org/10.1175/JCLI-D-14-00754.1)
-            and 
-            'deleted profile' OR Adrian Tompkins [email protected], posted on November 8, 2016 at
-            https://www.researchgate.net/post/Does-anyone-know-about-bias-correction-and-quantile-mapping-in-PYTHON
         '''
 
         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 == '+':
-            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))
-            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)
-
-            if kwargs.get('detrended', False):
-                '''detrended => shift mean of $T_{sim,p}$ to range of $T_{sim,h}$ to adjust extremes'''
-                for month, idxs in res.groupby('time.month').groups.items():
-                    m_simh, m_simp = [], []
-                    for idx in idxs:
-                        m_simh.append(simh[idx])
-                        m_simp.append(simp[idx])
-
-                    m_simh_mean = np.nanmean(m_simh)
-                    m_simp_mean = np.nanmean(m_simp)
-
-                    epsilon = np.interp((m_simp - m_simp_mean) + m_simh_mean, xbins, cdf_simh)
-                    X = (cls.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]
-            else:
-                epsilon = np.interp(simp, xbins, cdf_simh)                                  # Eq. 1
-                res.values = cls.get_inverse_of_cdf(cdf_obs, epsilon, xbins)                # Eq. 1
+        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))
+        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)
+
+        if kwargs.get('detrended', False) or kind in [ '*', 'mult' ]:
+            '''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 = [], []
+                for idx in idxs:
+                    m_simh.append(simh[idx])
+                    m_simp.append(simp[idx])
+                
+                m_simh = np.array(m_simh)
+                m_simp = np.array(m_simp)
+                m_simh_mean = np.nanmean(m_simh)     
+                m_simp_mean = np.nanmean(m_simp) 
+                
+                if kind == '+':
+                    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: 
+                    epsilon = np.interp((m_simh_mean * m_simp) / m_simp_mean, xbins, cdf_simh, left=.0, right=999.0)    # Eq. 2
+                    X = np.interp(epsilon, cdf_obs, xbins, left=.0, right=999.9) * (m_simp_mean / m_simh_mean)          # Eq. 2
+                    #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 
+            epsilon = np.interp(simp, xbins, cdf_simh)                                 # Eq. 1            
+            res.values = cls.get_inverse_of_cdf(cdf_obs, epsilon, xbins)               # Eq. 1      
             return res
-
-        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
-            '''
-
-            obs, simh, simp = np.sort(obs), np.sort(simh), np.array(simp)
-            global_max = max(np.amax(obs), np.amax(simh))
-            wide = global_max / n_quantiles
-            xbins = np.arange(0.0, global_max + wide, wide)
-
-            pdf_obs, bins = np.histogram(obs, bins=xbins)
-            pdf_simh, bins = np.histogram(simh, bins=xbins)
-
-            cdf_obs = np.insert(np.cumsum(pdf_obs), 0, 0.0)
-            cdf_simh = np.insert(np.cumsum(pdf_simh), 0, 0.0)
-
-            epsilon = np.interp(simp, xbins, cdf_simh, left=0.0, right=999.0)
-            return np.interp(epsilon, cdf_obs, xbins, left=0.0, right=-999.0)
-
         else: raise ValueError('Not implemented!')
 
     # ? -----========= E M P I R I C A L - Q U A N T I L E - M A P P I N G =========------
@@ -554,9 +552,7 @@ def empirical_quantile_mapping(cls,
         #         n_quantiles = n_quantiles,
         #         extrapolate = extrapolate
         #     )
-        # else:
-        #     obs, simh, simp = np.array(obs), np.array(simh), np.array(simp)
-        #     return corrected
+        # else: pass
 
     # ? -----========= Q U A N T I L E - D E L T A - M A P P I N G =========------
     @classmethod
@@ -638,7 +634,7 @@ def quantile_delta_mapping(cls,
             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
-            xbins = np.arange(kwargs.get('global_min', 0.0), global_max + wide, wide)
+            xbins = np.arange(kwargs.get('global_min', .0), global_max + wide, wide)
 
             cdf_obs = cls.get_cdf(obs,xbins)
             cdf_simh = cls.get_cdf(simh,xbins)
@@ -672,14 +668,14 @@ def get_inverse_of_cdf(base_cdf, insert_cdf, xbins) -> np.array:
     @staticmethod
     def load_data(
         obs_fpath: str,
-        contr_fpath: str,
-        scen_fpath: str,
+        simh_fpath: str,
+        simp_fpath: str,
         use_cftime: bool=False,
         chunks=None
     ) -> (xr.core.dataarray.Dataset, xr.core.dataarray.Dataset, xr.core.dataarray.Dataset):
         '''Load and return loaded netcdf datasets'''
         obs = xr.open_dataset(obs_fpath, use_cftime=use_cftime, chunks=chunks)
-        simh = xr.open_dataset(contr_fpath, use_cftime=use_cftime, chunks=chunks)
-        simp = xr.open_dataset(scen_fpath, use_cftime=use_cftime, chunks=chunks)
+        simh = xr.open_dataset(simh_fpath, use_cftime=use_cftime, chunks=chunks)
+        simp = xr.open_dataset(simp_fpath, use_cftime=use_cftime, chunks=chunks)
 
         return obs, simh, simp

examples/do_bias_correction.py

@@ -16,7 +16,7 @@
 # * ----- L O G G I N G -----
 formatter = logging.Formatter(
     fmt='%(asctime)s %(module)s,line: %(lineno)d %(levelname)8s | %(message)s',
-    datefmt='%Y-%m-%d %H:%M:%S' 
+    datefmt='%Y-%m-%d %H:%M:%S'
 )
 
 log = logging.getLogger()
@@ -37,7 +37,6 @@
 
 parser.add_argument('-g', '--group', dest='group', type=str, default=None, help='Value grouping, default: time, (options: time.month, time.dayofyear, time.year')
 parser.add_argument('-k', '--kind', dest='kind', type=str, default='+', help='+ or *, default: +')
-parser.add_argument('-w', '--window', dest='window', type=int, default=1, help='Window of grouping')
 parser.add_argument('-n', '--nquantiles', dest='n_quantiles', type=int, default=100, help='Nr. of Quantiles to use')
 
 parser.add_argument('-p', '--processes', dest='p', type=int, default=1, help='Multiprocessing with n processes, default: 1')
@@ -52,7 +51,6 @@
 unit = params['unit']
 group = params['group']
 kind = params['kind']
-window = params['window']
 n_quantiles = params['n_quantiles']
 n_jobs = params['p']
 
@@ -75,9 +73,8 @@ def main() -> None:
     end_date: str = ds_simp['time'][-1].dt.strftime('%Y%m%d').values.ravel()[0]
 
     descr1, descr2 = '', ''
-    if method in cm.XCLIM_SDBA_METHODS or method in ['quantile_mapping', 'quantile_delta_mapping']:
+    if method in cm.DISTRIBUTION_METHODS:
         descr1 = f'_quantiles-{n_quantiles}'
-        descr2 = f'_window-{window}'
 
     # ----- Adjustment -----
     log.info(f'Starting {method} adjustment')
@@ -89,7 +86,6 @@ def main() -> None:
         n_quantiles = n_quantiles,
         kind = kind,
         group = group,
-        window = window,
         n_jobs = n_jobs
     )
     log.info('Saving now')
@@ -103,4 +99,4 @@ def main() -> None:
     main()
 
 
-# * ----- ----- E O F ----- -----
+# * ----- ----- E O F ----- -----