-
Notifications
You must be signed in to change notification settings - Fork 2
/
sdm.py
170 lines (132 loc) · 6.78 KB
/
sdm.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
# BIAS CORRECTION
import numpy as np
import pandas as pd
from scipy.stats import gamma
from scipy.stats import norm
from scipy.signal import detrend
'''
Scaled distribution mapping for climate data
This is a excerpt from pyCAT and the method after Switanek et al. (2017) containing the functions to perform a relative and absolute bias correction on climate data.
(cc) c.jackisch@tu-braunschweig.de, July 2020
It is intended to be used on pandas time series at single locations/pixels.
Switanek, M. B., P. A. Troch, C. L. Castro, A. Leuprecht, H.-I. Chang, R. Mukherjee, and E. M. C. Demaria (2017), Scaled distribution mapping: a bias correction method that preserves raw climate model projected changes, Hydrol. Earth Syst. Sci., 21(6), 2649–2666, https://doi.org/10.5194/hess-21-2649-2017
'''
def relSDM(obs, mod, sce, cdf_threshold=0.9999999, lower_limit=0.1):
'''relative scaled distribution mapping assuming a gamma distributed parameter (with lower limit zero)
rewritten from pyCAT for 1D data
obs :: observed variable time series
mod :: modelled variable for same time series as obs
sce :: to unbias modelled time series
cdf_threshold :: upper and lower threshold of CDF
lower_limit :: lower limit of data signal (values below will be masked!)
returns corrected timeseries
tested with pandas series.
'''
obs_r = obs[obs >= lower_limit]
mod_r = mod[mod >= lower_limit]
sce_r = sce[sce >= lower_limit]
obs_fr = 1. * len(obs_r) / len(obs)
mod_fr = 1. * len(mod_r) / len(mod)
sce_fr = 1. * len(sce_r) / len(sce)
sce_argsort = np.argsort(sce)
obs_gamma = gamma.fit(obs_r, floc=0)
mod_gamma = gamma.fit(mod_r, floc=0)
sce_gamma = gamma.fit(sce_r, floc=0)
obs_cdf = gamma.cdf(np.sort(obs_r), *obs_gamma)
mod_cdf = gamma.cdf(np.sort(mod_r), *mod_gamma)
obs_cdf[obs_cdf > cdf_threshold] = cdf_threshold
mod_cdf[mod_cdf > cdf_threshold] = cdf_threshold
expected_sce_raindays = min(int(np.round(len(sce) * obs_fr * sce_fr / mod_fr)), len(sce))
sce_cdf = gamma.cdf(np.sort(sce_r), *sce_gamma)
sce_cdf[sce_cdf > cdf_threshold] = cdf_threshold
# interpolate cdf-values for obs and mod to the length of the scenario
obs_cdf_intpol = np.interp(np.linspace(1, len(obs_r), len(sce_r)), np.linspace(1, len(obs_r), len(obs_r)), obs_cdf)
mod_cdf_intpol = np.interp(np.linspace(1, len(mod_r), len(sce_r)), np.linspace(1, len(mod_r), len(mod_r)), mod_cdf)
# adapt the observation cdfs
obs_inverse = 1. / (1 - obs_cdf_intpol)
mod_inverse = 1. / (1 - mod_cdf_intpol)
sce_inverse = 1. / (1 - sce_cdf)
adapted_cdf = 1 - 1. / (obs_inverse * sce_inverse / mod_inverse)
adapted_cdf[adapted_cdf < 0.] = 0.
# correct by adapted observation cdf-values
xvals = gamma.ppf(np.sort(adapted_cdf), *obs_gamma) * gamma.ppf(sce_cdf, *sce_gamma) / gamma.ppf(sce_cdf,
*mod_gamma)
# interpolate to the expected length of future raindays
correction = np.zeros(len(sce))
if len(sce_r) > expected_sce_raindays:
xvals = np.interp(np.linspace(1, len(sce_r), expected_sce_raindays), np.linspace(1, len(sce_r), len(sce_r)),
xvals)
else:
xvals = np.hstack((np.zeros(expected_sce_raindays - len(sce_r)), xvals))
correction[sce_argsort[-expected_sce_raindays:]] = xvals
return pd.Series(correction, index=sce.index)
def absSDM(obs, mod, sce, cdf_threshold=0.9999999):
'''absolute scaled distribution mapping assuming a normal distributed parameter
rewritten from pyCAT for 1D data
obs :: observed variable time series
mod :: modelled variable for same time series as obs
sce :: to unbias modelled time series
cdf_threshold :: upper and lower threshold of CDF
returns corrected timeseries
tested with pandas series.
'''
obs_len = len(obs)
mod_len = len(mod)
sce_len = len(sce)
obs_mean = np.mean(obs)
mod_mean = np.mean(mod)
smean = np.mean(sce)
odetrend = detrend(obs)
mdetrend = detrend(mod)
sdetrend = detrend(sce)
obs_norm = norm.fit(odetrend)
mod_norm = norm.fit(mdetrend)
sce_norm = norm.fit(sdetrend)
sce_diff = sce - sdetrend
sce_argsort = np.argsort(sdetrend)
obs_cdf = norm.cdf(np.sort(odetrend), *obs_norm)
mod_cdf = norm.cdf(np.sort(mdetrend), *mod_norm)
sce_cdf = norm.cdf(np.sort(sdetrend), *sce_norm)
obs_cdf = np.maximum(np.minimum(obs_cdf, cdf_threshold), 1 - cdf_threshold)
mod_cdf = np.maximum(np.minimum(mod_cdf, cdf_threshold), 1 - cdf_threshold)
sce_cdf = np.maximum(np.minimum(sce_cdf, cdf_threshold), 1 - cdf_threshold)
# interpolate cdf-values for obs and mod to the length of the scenario
obs_cdf_intpol = np.interp(np.linspace(1, obs_len, sce_len), np.linspace(1, obs_len, obs_len), obs_cdf)
mod_cdf_intpol = np.interp(np.linspace(1, mod_len, sce_len), np.linspace(1, mod_len, mod_len), mod_cdf)
# adapt the observation cdfs
# split the tails of the cdfs around the center
obs_cdf_shift = obs_cdf_intpol - .5
mod_cdf_shift = mod_cdf_intpol - .5
sce_cdf_shift = sce_cdf - .5
obs_inverse = 1. / (.5 - np.abs(obs_cdf_shift))
mod_inverse = 1. / (.5 - np.abs(mod_cdf_shift))
sce_inverse = 1. / (.5 - np.abs(sce_cdf_shift))
adapted_cdf = np.sign(obs_cdf_shift) * (1. - 1. / (obs_inverse * sce_inverse / mod_inverse))
adapted_cdf[adapted_cdf < 0] += 1.
adapted_cdf = np.maximum(np.minimum(adapted_cdf, cdf_threshold), 1 - cdf_threshold)
xvals = norm.ppf(np.sort(adapted_cdf), *obs_norm) \
+ obs_norm[-1] / mod_norm[-1] \
* (norm.ppf(sce_cdf, *sce_norm) - norm.ppf(sce_cdf, *mod_norm))
xvals -= xvals.mean()
xvals += obs_mean + (smean - mod_mean)
correction = np.zeros(sce_len)
correction[sce_argsort] = xvals
correction += sce_diff - smean
return correction
def SDM(obs, mod, sce, meth='rel', cdf_threshold=0.9999999, lower_limit=0.1):
'''scaled distribution mapping - wrapper to relative and absolute bias correction functions
rewritten from pyCAT for 1D data
obs :: observed variable time series
mod :: modelled variable for same time series as obs
sce :: to unbias modelled time series
meth :: 'rel' for relative SDM, else absolute SDM will be performed
cdf_threshold :: upper and lower threshold of CDF
lower_limit :: lower limit of data signal (values below will be masked when meth != 'rel')
The original authors suggest to use the absolute SDM for air temperature and the relative SDM for precipitation and radiation series.
returns corrected timeseries
tested with pandas series.
'''
if meth == 'rel':
return relSDM(obs, mod, sce, cdf_threshold, lower_limit)
else:
return absSDM(obs, mod, sce, cdf_threshold)