Skip to content

Commit b4e4260

Browse files
adjusted and improved the examples; adjusted the readme
1 parent 75f2f22 commit b4e4260

File tree

7 files changed

+890
-288
lines changed

7 files changed

+890
-288
lines changed

README.md

Lines changed: 34 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -131,36 +131,57 @@ Notes:
131131

132132
Notebook with different methods and plots: `/examples/examples.ipynb`
133133

134-
Example script for adjusting climate data: `/examples/do_bias_correction.py`
134+
There is also an exmple script (`/examples/biasadjust.py`) that can be used to apply the available bias correction methods
135+
on 1- and 3-dimensional data sets (see `/examples/input_data/*.nc`).
136+
137+
Help:
135138

136139
```bash
137-
python3 do_bias_correction.py \
138-
--obs input_data/observations.nc \
140+
╰─ python3 biasadjust.py --help
141+
```
142+
143+
(1.) Example - Quantile Mapping bias correction on the provided example data:
144+
145+
```bash
146+
╰─ python3 biasadjust.py \
147+
--ref input_data/observations.nc \
139148
--contr input_data/control.nc \
140149
--scen input_data/scenario.nc \
141-
--method linear_scaling \
142-
--variable tas \
143-
--unit '°C' \
144-
--group 'time.month' \
145-
--kind +
150+
--kind "+" \
151+
--variable "tas" \
152+
-method quantile_mapping
146153
```
147154

148-
- Linear and variance, as well as delta change method require `--group time.month` as argument.
149-
- Adjustment methods that apply changes in distributional biases (QM, QDM, DQM, ...) require the `--nquantiles` argument set to some integer.
155+
(2.) Example - Linear Scaling bias correction on the provided example data:
156+
157+
```bash
158+
╰─ python3 biasadjust.py \
159+
--ref input_data/observations.nc \
160+
--contr input_data/control.nc \
161+
--scen input_data/scenario.nc \
162+
--kind "+" \
163+
--variable "tas" \
164+
--group "time.month" \
165+
-method linear_scaling
166+
```
167+
168+
Notes:
169+
150170
- Data sets must have the same spatial resolutions.
171+
- This script is far away from perfect - so please look at it, as a starting point. (:
151172

152173
---
153174

154175
<a name="notes"></a>
155176

156177
## 5. Notes
157178

158-
- 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).
159-
- Formulas and references can be found in the implementations of the corresponding functions.
179+
- 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).
180+
- 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/).
160181

161182
### Space for improvements:
162183

163-
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).
184+
- 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).
164185

165186
---
166187

cmethods/__init__.py

Lines changed: 26 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -141,6 +141,10 @@ def adjust_3d(
141141
"""
142142
Function to apply a bias correction method on 3-dimensional climate data.
143143
144+
*It is very important to pass ``group="time.month`` for scaling-based
145+
techniques if the correction should be performed as described in the
146+
referenced articles!* It is wanted to be the default.
147+
144148
:param method: The bias correction method to use
145149
:type method: str
146150
:param obs: The reference data set of the control period
@@ -214,7 +218,7 @@ def adjust_3d(
214218
"""
215219
obs = obs.transpose("lat", "lon", "time")
216220
simh = simh.transpose("lat", "lon", "time")
217-
simp = simp.transpose("lat", "lon", "time").load()
221+
simp = simp.transpose("lat", "lon", "time").load() # if not loaded ""
218222

219223
result = simp.copy(deep=True)
220224
len_lat, len_lon = len(simp.lat), len(simp.lon)
@@ -247,7 +251,7 @@ def adjust_3d(
247251
"simp": simp[lat],
248252
"group": group,
249253
"kind": kind,
250-
"n_quaniles": n_quantiles,
254+
"n_quantiles": n_quantiles,
251255
"kwargs": kwargs,
252256
}
253257
for lat in range(len_lat)
@@ -356,7 +360,7 @@ def linear_scaling(
356360
group: Union[str, None] = "time.month",
357361
kind: str = "+",
358362
**kwargs,
359-
) -> xr.core.dataarray.DataArray:
363+
) -> np.array:
360364
r"""
361365
The Linear Scaling bias correction technique can be applied on interval- and
362366
ratio-based climate variables to minimize deviations in the mean values
@@ -413,7 +417,7 @@ def linear_scaling(
413417
:type kind: str, optional
414418
:raises NotImplementedError: If the kind is not in (``+``, ``*``, ``add``, ``mult``)
415419
:return: The bias-corrected time series
416-
:rtype: xr.core.dataarray.DataArray
420+
:rtype: np.array
417421
418422
.. code-block:: python
419423
:linenos:
@@ -468,7 +472,7 @@ def variance_scaling(
468472
group: Union[str, None] = "time.month",
469473
kind: str = "+",
470474
**kwargs,
471-
) -> xr.core.dataarray.DataArray:
475+
) -> np.array:
472476
r"""
473477
The Variance Scaling bias correction technique can be applied on interval-based
474478
climate variables to minimize deviations in the mean and variance
@@ -533,7 +537,7 @@ def variance_scaling(
533537
:type kind: str, optional
534538
:raises NotImplementedError: If the kind is not in (``+``, ``add``)
535539
:return: The bias-corrected time series
536-
:rtype: xr.core.dataarray.DataArray
540+
:rtype: np.array
537541
538542
.. code-block:: python
539543
:linenos:
@@ -573,7 +577,7 @@ def variance_scaling(
573577
VS_1_simp = LS_simp - np.nanmean(LS_simp) # Eq. 4
574578

575579
adj_scaling_factor = cls.get_adjusted_scaling_factor(
576-
np.std(obs) / np.std(VS_1_simh),
580+
np.std(np.array(obs)) / np.std(VS_1_simh),
577581
kwargs.get("max_scaling_factor", cls.MAX_SCALING_FACTOR),
578582
)
579583

@@ -594,7 +598,7 @@ def delta_method(
594598
group: Union[str, None] = "time.month",
595599
kind: str = "+",
596600
**kwargs,
597-
) -> xr.core.dataarray.DataArray:
601+
) -> np.array:
598602
r"""
599603
The Delta Method bias correction technique can be applied on interval- and
600604
ratio-based climate variables to minimize deviations in the mean values
@@ -652,7 +656,7 @@ def delta_method(
652656
:type kind: str, optional
653657
:raises NotImplementedError: If the kind is not in (``+``, ``*``, ``add``, ``mult``)
654658
:return: The bias-corrected time series
655-
:rtype: xr.core.dataarray.DataArray
659+
:rtype: np.array
656660
657661
.. code-block:: python
658662
:linenos:
@@ -786,8 +790,7 @@ def quantile_mapping(
786790
... n_quantiles=250
787791
... )
788792
"""
789-
res = simp.copy(deep=True)
790-
obs, simh, simp = np.array(obs), np.array(simh), np.array(simp)
793+
obs, simh, simp_ = np.array(obs), np.array(simh), np.array(simp)
791794

792795
global_max = max(np.amax(obs), np.amax(simh))
793796
global_min = min(np.amin(obs), np.amin(simh))
@@ -799,7 +802,9 @@ def quantile_mapping(
799802

800803
if detrended:
801804
# detrended => shift mean of $X_{sim,p}$ to range of $X_{sim,h}$ to adjust extremes
802-
for _, idxs in res.groupby("time.month").groups.items():
805+
res = np.zeros(len(simp.values))
806+
for _, idxs in simp.groupby("time.month").groups.items():
807+
# detrended by long-term month
803808
m_simh, m_simp = [], []
804809
for idx in idxs:
805810
m_simh.append(simh[idx])
@@ -832,24 +837,22 @@ def quantile_mapping(
832837
m_simp_mean / m_simh_mean
833838
) # Eq. 2
834839
for i, idx in enumerate(idxs):
835-
res.values[idx] = X[i]
840+
res[idx] = X[i]
836841
return res
837842

838843
if kind in cls.ADDITIVE:
839-
epsilon = np.interp(simp, xbins, cdf_simh) # Eq. 1
840-
res.values = cls.get_inverse_of_cdf(cdf_obs, epsilon, xbins) # Eq. 1
841-
return res
844+
epsilon = np.interp(simp_, xbins, cdf_simh) # Eq. 1
845+
return cls.get_inverse_of_cdf(cdf_obs, epsilon, xbins) # Eq. 1
842846

843847
if kind in cls.MULTIPLICATIVE:
844848
epsilon = np.interp( # Eq. 2
845-
simp,
849+
simp_,
846850
xbins,
847851
cdf_simh,
848852
left=kwargs.get("val_min", 0.0),
849853
right=kwargs.get("val_max", None),
850854
)
851-
res.values = cls.get_inverse_of_cdf(cdf_obs, epsilon, xbins) # Eq. 2
852-
return res
855+
return cls.get_inverse_of_cdf(cdf_obs, epsilon, xbins) # Eq. 2
853856

854857
raise NotImplementedError(
855858
f"{kind} for quantile_mapping is not available. Use '+' or '*' instead."
@@ -992,7 +995,7 @@ def quantile_delta_mapping(
992995
:type kind: str, optional
993996
:raises NotImplementedError: If the kind is not in (``+``, ``*``, ``add``, ``mult``)
994997
:return: The bias-corrected time series
995-
:rtype: xr.core.dataarray.DataArray
998+
:rtype: np.array
996999
9971000
.. code-block:: python
9981001
:linenos:
@@ -1016,7 +1019,7 @@ def quantile_delta_mapping(
10161019
... )
10171020
"""
10181021
if kind in cls.ADDITIVE:
1019-
res = simp.copy(deep=True)
1022+
# res = simp.copy(deep=True)
10201023
obs, simh, simp = (
10211024
np.array(obs),
10221025
np.array(simh),
@@ -1035,11 +1038,9 @@ def quantile_delta_mapping(
10351038
epsilon = np.interp(simp, xbins, cdf_simp) # Eq. 1.1
10361039
QDM1 = cls.get_inverse_of_cdf(cdf_obs, epsilon, xbins) # Eq. 1.2
10371040
delta = simp - cls.get_inverse_of_cdf(cdf_simh, epsilon, xbins) # Eq. 1.3
1038-
res.values = QDM1 + delta # Eq. 1.4
1039-
return res
1041+
return QDM1 + delta # Eq. 1.4
10401042

10411043
if kind in cls.MULTIPLICATIVE:
1042-
res = simp.copy(deep=True)
10431044
obs, simh, simp = np.array(obs), np.array(simh), np.array(simp)
10441045
global_max = kwargs.get("global_max", max(np.amax(obs), np.amax(simh)))
10451046
wide = global_max / n_quantiles
@@ -1058,13 +1059,11 @@ def quantile_delta_mapping(
10581059
) # Eq. 2.3
10591060
delta[np.isnan(delta)] = 0
10601061

1061-
res.values = QDM1 * delta # Eq. 2.4
1062-
return res
1062+
return QDM1 * delta # Eq. 2.4
10631063
raise NotImplementedError(
10641064
f"{kind} not available for quantile_delta_mapping. Use '+' or '*' instead."
10651065
)
10661066

1067-
# * -----========= G E N E R A L =========------
10681067
@staticmethod
10691068
def get_pdf(x: Union[list, np.array], xbins: Union[list, np.array]) -> np.array:
10701069
r"""

0 commit comments

Comments
 (0)