-
Notifications
You must be signed in to change notification settings - Fork 2
/
DAA_EnKF_batchOSS_phenoShift.py
375 lines (323 loc) · 16.9 KB
/
DAA_EnKF_batchOSS_phenoShift.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
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
# -*- coding: utf-8 -*-
"""
Created on Wed Jul 21 15:17:52 2021
@author: Qi Yang
OSS Experiment for only assimilating Canopy Cover (CC)
Evaluate the model performance at a wrong guess of planting date
"""
import numpy as np
import matlab.engine
from KFs import EnsembleKalmanFilter as EnKF
from DAA_UncertainPara import uncertain_para
import matplotlib.pyplot as plt
import matplotlib
import os
matplotlib.rcParams['font.family'] = 'Times New Roman'
matplotlib.rcParams['figure.dpi'] = 300
class Aquacrop_env():
def __init__(self, inputpath, start_date, ensemble_n, init_para, state_case = 1, initUpdate = False):
"""
Parameters
----------
inputpath : str
output : str
start_date : list
ensemble_n : int
init_para : list >> the list of n (sample number) * N (parameter number) parameters
"""
self.eng = matlab.engine.start_matlab()
self.eng.addpath('AquaCropOS_v60a')
self.inputpath = inputpath
self.start_date = start_date
self.ensemble_n = ensemble_n
self.init_para = init_para
if len(init_para) != ensemble_n:
raise ValueError('ensemble samples number error!')
self.init_state = []
self.state_case = state_case
self.updateNextStep = initUpdate
if self.state_case == 1:
self.stateList = ['CC']
elif self.state_case == 2:
self.stateList = ['CC', 'biomass']
elif self.state_case == 3:
self.stateList = ['CC', 'biomass', 'Para_CDC', 'Para_CGC','Para_CCx']
elif self.state_case == 4:
self.stateList = ['CC', 'biomass', 'Para_CDC', 'Para_CGC','Para_CCx', 'Para_WP','Para_Zmax','Para_p_up2','Para_p_up4']
elif self.state_case == 5:
self.stateList = ['CC', 'biomass', 'Para_CDC', 'Para_CGC','Para_CCx', 'Para_WP','Para_Zmax','Para_p_up2','Para_p_up4'
,'GDDcum','Para_HIstart','Para_YldForm', 'Para_Senescence']
elif self.state_case == 6:
self.stateList = ['CC', 'biomass', 'Para_CDC', 'Para_CGC','Para_CCx', 'Para_WP','Para_Zmax','Para_p_up2','Para_p_up4'
, 'Para_Tbase','Para_HIstart','Para_YldForm', 'Para_Senescence']
elif self.state_case == 7: # update the planting date and other revelated paras based on GDD
self.stateList = ['CC', 'biomass', 'Para_CDC', 'Para_CGC','Para_CCx', 'Para_WP','Para_Zmax','Para_p_up2','Para_p_up4'
, 'GDDcum','Para_HIstart','Para_YldForm', 'Para_Senescence']
elif self.state_case == 8:
self.stateList = ['CC', 'biomass', 'Para_CDC', 'Para_CGC','Para_CCx','Para_WP','Para_Zmax','Para_p_up2','Para_p_up4',
'GDDcum','Para_Tbase','Para_HIstart','Para_YldForm', 'Para_Senescence']
self.reset()
def reset(self):
print ('Initializing the ensemble models...')
for n in range(self.ensemble_n):
# engine start and initialaze
if self.init_para != None:
self.eng.DA_pyaqua_init(inputpath, start_date,
matlab.double(self.init_para[n]), int(n + 1), nargout= 0) # inside matlab.double should be a list
else:
self.eng.DA_pyaqua_init(inputpath, start_date,
[], int(n + 1), nargout= 0)
def steprun(self, state_in, dt, sample_n):
out = self.eng.step_run_DA_OSS(matlab.double(list(state_in)), int(dt), int(sample_n + 1), int(self.state_case), self.updateNextStep) # sample_n start from 1
state_out = [out[t] for t in self.stateList]
return state_out, out
def hx(x):
return np.array([x[0]])
def plotState(xs_enkf, color = 'r', DA_state_loc = 0, P_enkf = None, label = ' '):
y = np.asarray([x[DA_state_loc] for x in xs_enkf])
plt.plot(range(len(y)),y,color = color, label = label)
ps = [p[DA_state_loc,DA_state_loc]**.5 for p in P_enkf]
ps = np.asarray(ps)
plt.fill_between(range(len(y)), y + ps ,
y - ps,color = color,alpha=.25,linewidth=0.5,linestyle='--')
def plotState_extent(All_state, state_key, color = 'g', color2 = 'darkgreen', label = ' '):
y = []
std = []
Yall=[]
for state in All_state:
tmp = []
for s in state:
tmp.append(s[state_key])
Yall.append(tmp)
std.append(np.std(tmp,ddof = 1)) # ddof = 1 mean no-bias std, divide (n-1)
y.append(np.mean(tmp))
y = np.asarray(y)
Yall = np.array(Yall)
for i in range(Yall.shape[1]):
plt.plot(range(len(Yall[:,i])),Yall[:,i],color = color,alpha=.10,linewidth=0.5)
plt.plot(range(len(y)),y,color = color2, label = label)
# plt.plot(range(len(y)),y,color = color, label = label)
# plt.fill_between(range(len(y)), y + std ,
# y - std,color = color,alpha=.25,linewidth=0.5,linestyle='--')
def DA_plot(state_key, y_trueth, ylabel, title, plotPara = False):
# fig = plt.figure()
fig, ax = plt.subplots(1, 1,figsize = (6,5))
plotState_extent(All_state, state_key = state_key, color = 'b', color2 = 'darkblue',label = 'Open-loop')
plotState_extent(All_state2, state_key = state_key, color = 'g', color2 = 'darkgreen', label = 'EnKF')
if plotPara:
plt.axhline(y_trueth, color = 'k', linestyle = '--', label = 'Reference')
else:
plt.plot(range(len(y_trueth)),y_trueth, color = 'k', linestyle = '--', label = 'Reference')
plt.xlim(0,100)
plt.legend()
plt.xlabel('Day after plant',fontsize=16)
plt.ylabel(ylabel,fontsize=16)
# plt.title(title)
plt.text(0.45, 0.9, title, transform=ax.transAxes,fontsize=16)
return fig
def DA_plot_state(DA_state_loc, y_trueth, ylabel, title):
fig = plt.figure()
plotState(xs_enkf = xs_enkf, color = 'r', DA_state_loc = DA_state_loc, P_enkf = P_enkf, label = 'Open-loop')
plotState(xs_enkf = xs_enkf2, color = 'g', DA_state_loc = DA_state_loc, P_enkf = P_enkf2, label = 'EnKF')
# plotState_extent(All_state2, state_key = 'CC', color = 'b', DA_state_loc = 0, P_enkf = P_enkf2, label = 'EnKF-modelstate')
plt.plot(range(len(y_trueth)),y_trueth, color = 'k', label = 'Reference')
plt.xlim(0,100)
plt.legend()
plt.xlabel('Day after plant')
plt.ylabel(ylabel)
plt.title(title)
return fig
if __name__ == '__main__':
## model settings
inputpath = './AquaCropOS_v60a/Input/rice_NN/'
# the planting date was moved ahead by 10 days
aheadDay = 10
start_date = ['2039','7','%d'%(15-aheadDay)]
ensemble_n = 30
np.random.seed(0)
outPath = 'OSSresult_1obs_phenoShift'#'OSSresult_1obs'
saveResult = False
forceOpenloop = False
OSS_path='OSSE_truth'
## Ensemble samples sampling
UP = uncertain_para()
u_para_default = UP.u_para_default
CV = [0.5, 0.05, 0.05, 0.08, 0.1,
0.1, 0.1, 0.05, 0.05, 0.1,
0.05, 0.15, 0.1, 0.1]
std2_0 = list((np.asarray(u_para_default) * np.asarray(CV))**2)
P_u_para = np.diag(std2_0)
init_paras = list(np.random.multivariate_normal(mean=u_para_default, cov=P_u_para, size=ensemble_n))
for i,init_para in enumerate(init_paras):
extent_para = UP.replace_uncertain(init_para) # padding the uncertain parameters to full size
valid_para = UP.validation_check(extent_para) # correct the init parameters
init_paras[i] = valid_para
## batch Run EnKF
# state_case_list = [1,2,3,4,5,6]
state_case_list = [8]
# load the "ground truth" of OSS experiment which was generated by one random run
z_origin = list(np.zeros(aheadDay))+list(np.loadtxt('%s/sim_z.in'%OSS_path))
z_para = np.loadtxt('%s/sim_para.in'%OSS_path)
for case_n, state_case in enumerate(state_case_list):
print ('Runing case %d ...'%state_case)
## initialized the ensembled Aquacrop models
Aquacrop = Aquacrop_env(inputpath, start_date, ensemble_n, init_para = init_paras, state_case = state_case, initUpdate = False)
## EnKF settings
x_paraLoc = None # If the first step needs to update and discrete the states, should give the paralocation and set initUpdate = True
R = np.diag([0.03**2]) # measurement cov matrix
P0 = np.diag([0.])
if case_n == 0:
if (not os.path.exists('openLoopResultShift.npy')) | forceOpenloop:
# ENKF - open loop
DAA_enkf = EnKF(x=np.zeros(len(Aquacrop.stateList)), P=P0, dim_z=1, N=ensemble_n, hx = hx, fx = Aquacrop.steprun, x_paraLoc = x_paraLoc)
DAA_enkf.R = R
dt = 1
xs_enkf = []
sigmas_enkf = []
P_enkf = []
K_enkf = []
All_state = []
count = 0
while not DAA_enkf.allModelDone:
# predict
DAA_enkf.predict(dt = dt)
xs_enkf.append(DAA_enkf.x.copy())
P_enkf.append(DAA_enkf.P.copy())
K_enkf.append(DAA_enkf.K.copy())
sigmas_enkf.append(DAA_enkf.sigmas.copy())
All_state.append(DAA_enkf.allState.copy())
count += 1
if count%10 == 0:
print ('Simulating %d days...'%count)
Aquacrop.reset()
np.save('openLoopResultShift.npy',[xs_enkf,P_enkf,K_enkf,sigmas_enkf,All_state])
else:
xs_enkf,P_enkf,K_enkf,sigmas_enkf,All_state = np.load('openLoopResultShift.npy',allow_pickle=True)
# ENKF - assimilation
DAA_enkf = EnKF(x=np.zeros(len(Aquacrop.stateList)), P=P0, dim_z=1, N=ensemble_n, hx = hx, fx = Aquacrop.steprun, x_paraLoc = x_paraLoc)
DAA_enkf.R = R.copy()
dt = 1
xs_enkf2 = []
sigmas_enkf2 = []
P_enkf2 = []
K_enkf2 = []
All_state2 = []
z_DAP = list(np.arange(1,len(z_origin),5))
z = [t for i,t in enumerate(z_origin) if (i+1) in z_DAP]
DAP = 1
ass_n = 0
while not DAA_enkf.allModelDone:
# predict
DAA_enkf.predict(dt = dt)
Aquacrop.updateNextStep = False
# update
if DAP in z_DAP:
DAA_enkf.update(np.asarray([z[ass_n]]))
Aquacrop.updateNextStep = True
ass_n += 1
if DAP%10 == 0:
print ('Simulating %d days...'%DAP)
# record
xs_enkf2.append(DAA_enkf.x.copy())
P_enkf2.append(DAA_enkf.P.copy())
K_enkf2.append(DAA_enkf.K.copy())
sigmas_enkf2.append(DAA_enkf.sigmas.copy())
All_state2.append(DAA_enkf.allState.copy())
DAP += 1
print ('Simulation of case %d is finished...'%state_case)
if saveResult:
if not os.path.exists(outPath):
os.mkdir(outPath)
np.save('%s/case%d.npy'%(outPath,state_case),[All_state,All_state2])
# plot CC
# fig = DA_plot_state(DA_state_loc = 0, y_trueth = z_origin, ylabel = 'Canoy cover %', title = 'Case %d'%state_case)
fig = DA_plot(state_key = 'CC', y_trueth = z_origin, ylabel = 'Canoy cover %', title = 'Case %d'%state_case)
if saveResult:
fig.savefig(outPath + '/CC_case%d.jpg'%state_case, bbox_inches='tight')
# plot yield
y_trueth = list(np.zeros(aheadDay))+list(np.loadtxt('%s/sim_yield.in'%OSS_path))
fig = DA_plot(state_key = 'yield',y_trueth = y_trueth, ylabel = 'Yield t/ha', title = 'Case %d'%state_case)
plt.ylim(-0.2,8)
if saveResult:
fig.savefig(outPath + '/Yield_case%d.jpg'%state_case, bbox_inches='tight')
# plot Biomass
y_trueth = list(np.zeros(aheadDay))+list(np.loadtxt('%s/sim_bio.in'%OSS_path))
fig = DA_plot(state_key = 'biomass',y_trueth = y_trueth, ylabel = 'Biomass g/m2', title = 'Case %d'%state_case)
plt.ylim(-50,2000)
if saveResult:
fig.savefig(outPath + '/Bio_case%d.jpg'%state_case, bbox_inches='tight')
if state_case == 6:
# plot Tbase
y_trueth = z_para[21]
fig = DA_plot(state_key = 'Para_Tbase',y_trueth = y_trueth, ylabel = 'Tbase C', plotPara = True, title = 'Case %d'%state_case)
plt.ylim(7,11)
if saveResult:
fig.savefig(outPath + '/Tbase_case%d.jpg'%state_case, bbox_inches='tight')
if state_case in [5,6]:
# plot phenology parameters
y_trueth = z_para[11]
fig = DA_plot(state_key = 'Para_HIstart',y_trueth = y_trueth, ylabel = 'HIstart', plotPara = True, title = 'Case %d'%state_case)
plt.ylim(1000,1300)
if saveResult:
fig.savefig(outPath + '/HIstart_case%d.jpg'%state_case, bbox_inches='tight')
y_trueth = z_para[14]
fig = DA_plot(state_key = 'Para_YldForm',y_trueth = y_trueth, ylabel = 'YldForm', plotPara = True, title = 'Case %d'%state_case)
plt.ylim(600,850)
if saveResult:
fig.savefig(outPath + '/YldForm_case%d.jpg'%state_case, bbox_inches='tight')
y_trueth = z_para[12]
fig = DA_plot(state_key = 'Para_Senescence',y_trueth = y_trueth, ylabel = 'Senescence', plotPara = True, title = 'Case %d'%state_case)
plt.ylim(1000,1400)
if saveResult:
fig.savefig(outPath + '/Sen_case%d.jpg'%state_case, bbox_inches='tight')
if state_case >= 3:
# plot CDC
y_trueth = z_para[20]
fig = DA_plot(state_key = 'Para_CDC',y_trueth = y_trueth, ylabel = 'CDC', plotPara = True, title = 'Case %d'%state_case)
plt.ylim(0.003,0.006)
if saveResult:
fig.savefig(outPath + '/CDC_case%d.jpg'%state_case, bbox_inches='tight')
# plot CGC
y_trueth = z_para[18]
fig = DA_plot(state_key = 'Para_CGC',y_trueth = y_trueth, ylabel = 'CGC', plotPara = True, title = 'Case %d'%state_case)
plt.ylim(0.003,0.009)
if saveResult:
fig.savefig(outPath + '/CGC_case%d.jpg'%state_case, bbox_inches='tight')
# plot CCx
y_trueth = z_para[19]
fig = DA_plot(state_key = 'Para_CCx',y_trueth = y_trueth, ylabel = 'CCx', plotPara = True, title = 'Case %d'%state_case)
plt.ylim(0.6,1)
if saveResult:
fig.savefig(outPath + '/CCx_case%d.jpg'%state_case, bbox_inches='tight')
if state_case >= 4:
# plot WP
y_trueth = z_para[32]
fig = DA_plot(state_key = 'Para_WP',y_trueth = y_trueth, ylabel = 'WP', plotPara = True, title = 'Case %d'%state_case)
plt.ylim(14,21)
if saveResult:
fig.savefig(outPath + '/WP_case%d.jpg'%state_case, bbox_inches='tight')
# plot Zmax
y_trueth = z_para[24]
fig = DA_plot(state_key = 'Para_Zmax',y_trueth = y_trueth, ylabel = 'Zmax', plotPara = True, title = 'Case %d'%state_case)
plt.ylim(0.3, 0.6)
if saveResult:
fig.savefig(outPath + '/Zmax_case%d.jpg'%state_case, bbox_inches='tight')
# plot p_up2
y_trueth = z_para[39]
fig = DA_plot(state_key = 'Para_p_up2',y_trueth = y_trueth, ylabel = 'p_up2', plotPara = True, title = 'Case %d'%state_case)
plt.ylim(0.3, 0.6)
if saveResult:
fig.savefig(outPath + '/p_up2_case%d.jpg'%state_case, bbox_inches='tight')
# plot p_up4
y_trueth = z_para[43]
fig = DA_plot(state_key = 'Para_p_up4',y_trueth = y_trueth, ylabel = 'p_up4', plotPara = True, title = 'Case %d'%state_case)
plt.ylim(0.4, 0.9)
if saveResult:
fig.savefig(outPath + '/p_up4_case%d.jpg'%state_case, bbox_inches='tight')
if state_case in [5,6]:
# plot GDD
y_trueth = list(np.zeros(aheadDay))+list(np.loadtxt('%s/sim_GDD.in'%OSS_path))
fig = DA_plot(state_key = 'GDDcum', y_trueth = y_trueth, ylabel = 'GDD g/m2', title = 'Case %d'%state_case)
plt.ylim(-50,2200)
if saveResult:
fig.savefig(outPath + '/GDD_case%d.jpg'%state_case, bbox_inches='tight')