Merge pull request #15 from btschwertfeger/adjust-examples

Adjusted and improved the examples

Author: btschwertfeger

README.md

@@ -1,4 +1,4 @@
-# Bias correction/adjustment procedures for climatic reasearch
+# python-cmethods
 
 <div align="center">
 
@@ -18,7 +18,7 @@
 
 </div>
 
-This Python module contains a collection of different scale- and distribution-based bias adjustment techniques for climatic research (see `/examples/examples.ipynb` for help).
+This Python module serves as a collection of different scale- and distribution-based bias correction techniques for climatic research
 
 The documentation is available at: [https://python-kraken-sdk.readthedocs.io/en/stable/](https://python-kraken-sdk.readthedocs.io/en/stable/)
 
@@ -67,7 +67,7 @@ In this way, for example, modeled data, which on average represent values that a
 
 ## 2. Available methods
 
-All methods except the `adjust_3d` function requires the application on one time series.
+All methods except the `adjust_3d` function requires that the input data sets only contain one dimension.
 
 | Function name            | Description                                                                                                                                                  |
 | ------------------------ | ------------------------------------------------------------------------------------------------------------------------------------------------------------ |
@@ -131,36 +131,57 @@ Notes:
 
 Notebook with different methods and plots: `/examples/examples.ipynb`
 
-Example script for adjusting climate data: `/examples/do_bias_correction.py`
+There is also an exmple script (`/examples/biasadjust.py`) that can be used to apply the available bias correction methods
+on 1- and 3-dimensional data sets (see `/examples/input_data/*.nc`).
+
+Help:
 
 ```bash
-python3 do_bias_correction.py         \
-    --obs input_data/observations.nc  \
+╰─ python3 biasadjust.py --help
+```
+
+(1.) Example - Quantile Mapping bias correction on the provided example data:
+
+```bash
+╰─ python3 biasadjust.py              \
+    --ref input_data/observations.nc  \
     --contr input_data/control.nc     \
     --scen input_data/scenario.nc     \
-    --method linear_scaling           \
-    --variable tas                    \
-    --unit '°C'                       \
-    --group 'time.month'              \
-    --kind +
+    --kind "+"                        \
+    --variable "tas"                  \
+    --method quantile_mapping
 ```
 
-- Linear and variance, as well as delta change method require `--group time.month` as argument.
-- Adjustment methods that apply changes in distributional biases (QM, QDM, DQM, ...) require the `--nquantiles` argument set to some integer.
+(2.) Example - Linear Scaling bias correction on the provided example data:
+
+```bash
+╰─ python3 biasadjust.py              \
+    --ref input_data/observations.nc  \
+    --contr input_data/control.nc     \
+    --scen input_data/scenario.nc     \
+    --kind "+"                        \
+    --variable "tas"                  \
+    --group "time.month"              \
+    --method linear_scaling
+```
+
+Notes:
+
 - Data sets must have the same spatial resolutions.
+- This script is far away from perfect - so please look at it, as a starting point. (:
 
 ---
 
 <a name="notes"></a>
 
 ## 5. Notes
 
-- Computation in Python takes some time, so this is only for demonstration. When adjusting large datasets, its best to use the C++ tool [BiasAdjustCXX](https://github.com/btschwertfeger/BiasAdjustCXX).
-- Formulas and references can be found in the implementations of the corresponding functions.
+- Computation in Python takes some time, so this is only for demonstration. When adjusting large datasets, its best to use the command-line tool [BiasAdjustCXX](https://github.com/btschwertfeger/BiasAdjustCXX).
+- Formulas and references can be found in the implementations of the corresponding functions, on the bottom of the README.md and in the [documentation](https://python-kraken-sdk.readthedocs.io/en/stable/).
 
 ### Space for improvements:
 
-Since the scaling methods implemented so far scale by default over the mean values of the respective months, unrealistic long-term mean values may occur at the month transitions. This can be prevented either by selecting `group='time.dayofyear'`. Alternatively, it is possible not to scale using long-term mean values, but using a 31-day interval, which takes the 31 surrounding values over all years as the basis for calculating the mean values. This is not yet implemented in this module, but is available in the C++ tool [BiasAdjustCXX](https://github.com/btschwertfeger/BiasAdjustCXX).
+- Since the scaling methods implemented so far scale by default over the mean values of the respective months, unrealistic long-term mean values may occur at the month transitions. This can be prevented either by selecting `group='time.dayofyear'`. Alternatively, it is possible not to scale using long-term mean values, but using a 31-day interval, which takes the 31 surrounding values over all years as the basis for calculating the mean values. This is not yet implemented, because even the computation for this takes so much time, that it is not worth implementing it in python - but this is available in [BiasAdjustCXX](https://github.com/btschwertfeger/BiasAdjustCXX).
 
 ---
 

cmethods/__init__.py

@@ -53,7 +53,7 @@ class CMethods:
     The following bias correction techniques are available:
         Scaling-based techniques:
             * Linear Scaling :func:`cmethods.CMethods.linear_scaling`
-            * Vairance Scaling :func:`cmethods.CMethods.variance_scaling`
+            * Variance Scaling :func:`cmethods.CMethods.variance_scaling`
             * Delta (change) Method :func:`cmethods.CMethods.delta_method`
 
         Distribution-based techniques:
@@ -137,10 +137,14 @@ def adjust_3d(
         group: Union[str, None] = None,
         n_jobs: int = 1,
         **kwargs,
-    ) -> xr.core.dataarray.Dataset:
+    ) -> xr.core.dataarray.DataArray:
         """
         Function to apply a bias correction method on 3-dimensional climate data.
 
+        *It is very important to pass ``group="time.month`` for scaling-based
+        techniques if the correction should be performed as described in the
+        referenced articles!* It is wanted to be the default.
+
         :param method: The bias correction method to use
         :type method: str
         :param obs: The reference data set of the control period
@@ -162,7 +166,7 @@ def adjust_3d(
         :type n_jobs: int, optional
         :raises UnknownMethodError: If the correction method is not implemented
         :return: The bias-corrected time series
-        :rtype: xr.core.dataarray.Dataset
+        :rtype: xr.core.dataarray.DataArray
 
         .. code-block:: python
             :linenos:
@@ -247,7 +251,7 @@ def adjust_3d(
                             "simp": simp[lat],
                             "group": group,
                             "kind": kind,
-                            "n_quaniles": n_quantiles,
+                            "n_quantiles": n_quantiles,
                             "kwargs": kwargs,
                         }
                         for lat in range(len_lat)
@@ -344,7 +348,7 @@ def grouped_correction(
             for i, index in enumerate(groups[month]):
                 result[index] = computed_result[i]
 
-        return result
+        return np.array(result)
 
     # ? -----========= L I N E A R - S C A L I N G =========------
     @classmethod
@@ -356,7 +360,7 @@ def linear_scaling(
         group: Union[str, None] = "time.month",
         kind: str = "+",
         **kwargs,
-    ) -> xr.core.dataarray.DataArray:
+    ) -> np.array:
         r"""
         The Linear Scaling bias correction technique can be applied on interval- and
         ratio-based climate variables to minimize deviations in the mean values
@@ -413,7 +417,7 @@ def linear_scaling(
         :type kind: str, optional
         :raises NotImplementedError: If the kind is not in (``+``, ``*``, ``add``, ``mult``)
         :return: The bias-corrected time series
-        :rtype: xr.core.dataarray.DataArray
+        :rtype: np.array
 
         .. code-block:: python
             :linenos:
@@ -468,7 +472,7 @@ def variance_scaling(
         group: Union[str, None] = "time.month",
         kind: str = "+",
         **kwargs,
-    ) -> xr.core.dataarray.DataArray:
+    ) -> np.array:
         r"""
         The Variance Scaling bias correction technique can be applied on interval-based
         climate variables to minimize deviations in the mean and variance
@@ -533,7 +537,7 @@ def variance_scaling(
         :type kind: str, optional
         :raises NotImplementedError: If the kind is not in (``+``, ``add``)
         :return: The bias-corrected time series
-        :rtype: xr.core.dataarray.DataArray
+        :rtype: np.array
 
         .. code-block:: python
             :linenos:
@@ -573,7 +577,7 @@ def variance_scaling(
             VS_1_simp = LS_simp - np.nanmean(LS_simp)  # Eq. 4
 
             adj_scaling_factor = cls.get_adjusted_scaling_factor(
-                np.std(obs) / np.std(VS_1_simh),
+                np.std(np.array(obs)) / np.std(VS_1_simh),
                 kwargs.get("max_scaling_factor", cls.MAX_SCALING_FACTOR),
             )
 
@@ -594,7 +598,7 @@ def delta_method(
         group: Union[str, None] = "time.month",
         kind: str = "+",
         **kwargs,
-    ) -> xr.core.dataarray.DataArray:
+    ) -> np.array:
         r"""
         The Delta Method bias correction technique can be applied on interval- and
         ratio-based climate variables to minimize deviations in the mean values
@@ -652,7 +656,7 @@ def delta_method(
         :type kind: str, optional
         :raises NotImplementedError: If the kind is not in (``+``, ``*``, ``add``, ``mult``)
         :return: The bias-corrected time series
-        :rtype: xr.core.dataarray.DataArray
+        :rtype: np.array
 
         .. code-block:: python
             :linenos:
@@ -707,7 +711,7 @@ def quantile_mapping(
         kind: str = "+",
         detrended: bool = False,
         **kwargs,
-    ) -> xr.core.dataarray.DataArray:
+    ) -> np.array:
         r"""
         The Quantile Mapping bias correction technique can be used to minimize distributional
         biases between modeled and observed time-series climate data. Its interval-independant
@@ -763,7 +767,7 @@ def quantile_mapping(
         :type detrended: bool, optional
         :raises NotImplementedError: If the kind is not in (``+``, ``*``, ``add``, ``mult``)
         :return: The bias-corrected time series
-        :rtype: xr.core.dataarray.DataArray
+        :rtype: np.array
 
         .. code-block:: python
             :linenos:
@@ -786,8 +790,7 @@ def quantile_mapping(
             ...     n_quantiles=250
             ... )
         """
-        res = simp.copy(deep=True)
-        obs, simh, simp = np.array(obs), np.array(simh), np.array(simp)
+        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))
@@ -799,7 +802,9 @@ def quantile_mapping(
 
         if detrended:
             # detrended => shift mean of $X_{sim,p}$ to range of $X_{sim,h}$ to adjust extremes
-            for _, idxs in res.groupby("time.month").groups.items():
+            res = np.zeros(len(simp.values))
+            for _, idxs in simp.groupby("time.month").groups.items():
+                # detrended by long-term month
                 m_simh, m_simp = [], []
                 for idx in idxs:
                     m_simh.append(simh[idx])
@@ -832,24 +837,22 @@ def quantile_mapping(
                         m_simp_mean / m_simh_mean
                     )  # Eq. 2
                 for i, idx in enumerate(idxs):
-                    res.values[idx] = X[i]
+                    res[idx] = X[i]
             return res
 
         if kind in cls.ADDITIVE:
-            epsilon = np.interp(simp, xbins, cdf_simh)  # Eq. 1
-            res.values = cls.get_inverse_of_cdf(cdf_obs, epsilon, xbins)  # Eq. 1
-            return res
+            epsilon = np.interp(simp_, xbins, cdf_simh)  # Eq. 1
+            return cls.get_inverse_of_cdf(cdf_obs, epsilon, xbins)  # Eq. 1
 
         if kind in cls.MULTIPLICATIVE:
             epsilon = np.interp(  # Eq. 2
-                simp,
+                simp_,
                 xbins,
                 cdf_simh,
                 left=kwargs.get("val_min", 0.0),
                 right=kwargs.get("val_max", None),
             )
-            res.values = cls.get_inverse_of_cdf(cdf_obs, epsilon, xbins)  # Eq. 2
-            return res
+            return cls.get_inverse_of_cdf(cdf_obs, epsilon, xbins)  # Eq. 2
 
         raise NotImplementedError(
             f"{kind} for quantile_mapping is not available. Use '+' or '*' instead."
@@ -900,7 +903,7 @@ def quantile_delta_mapping(
         n_quantiles: int,
         kind: str = "+",
         **kwargs,
-    ) -> xr.core.dataarray.DataArray:
+    ) -> np.array:
         r"""
         The Quantile Delta Mapping bias correction technique can be used to minimize distributional
         biases between modeled and observed time-series climate data. Its interval-independant
@@ -992,7 +995,7 @@ def quantile_delta_mapping(
         :type kind: str, optional
         :raises NotImplementedError: If the kind is not in (``+``, ``*``, ``add``, ``mult``)
         :return: The bias-corrected time series
-        :rtype: xr.core.dataarray.DataArray
+        :rtype: np.array
 
         .. code-block:: python
             :linenos:
@@ -1016,7 +1019,6 @@ def quantile_delta_mapping(
             ... )
         """
         if kind in cls.ADDITIVE:
-            res = simp.copy(deep=True)
             obs, simh, simp = (
                 np.array(obs),
                 np.array(simh),
@@ -1035,11 +1037,9 @@ def quantile_delta_mapping(
             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
-            res.values = QDM1 + delta  # Eq. 1.4
-            return res
+            return QDM1 + delta  # Eq. 1.4
 
         if kind in cls.MULTIPLICATIVE:
-            res = simp.copy(deep=True)
             obs, simh, simp = np.array(obs), np.array(simh), np.array(simp)
             global_max = kwargs.get("global_max", max(np.amax(obs), np.amax(simh)))
             wide = global_max / n_quantiles
@@ -1058,13 +1058,11 @@ def quantile_delta_mapping(
                 )  # Eq. 2.3
             delta[np.isnan(delta)] = 0
 
-            res.values = QDM1 * delta  # Eq. 2.4
-            return res
+            return QDM1 * delta  # Eq. 2.4
         raise NotImplementedError(
             f"{kind} not available for quantile_delta_mapping. Use '+' or '*' instead."
         )
 
-    # * -----========= G E N E R A L  =========------
     @staticmethod
     def get_pdf(x: Union[list, np.array], xbins: Union[list, np.array]) -> np.array:
         r"""

docs/src/introduction.rst

@@ -45,7 +45,7 @@ Available Methods
 The following bias correction techniques are available:
     Scaling-based techniques:
         * Linear Scaling :func:`cmethods.CMethods.linear_scaling`
-        * Vairance Scaling :func:`cmethods.CMethods.variance_scaling`
+        * Variance Scaling :func:`cmethods.CMethods.variance_scaling`
         * Delta (change) Method :func:`cmethods.CMethods.delta_method`
 
     Distribution-based techniques: