-
Notifications
You must be signed in to change notification settings - Fork 0
/
0_1_calRMSE_freerun_da.py
116 lines (105 loc) · 3.11 KB
/
0_1_calRMSE_freerun_da.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
# -*- coding: utf-8 -*-
"""
Created on Tue Jul 19 15:17:38 2022
@author: Zikang He
"""
import numpy as np
import matplotlib.pyplot as plt
from math import sqrt
xx = np.load('true_l.npy')
E = np.load('E.npz')
da =E['E']
E = np.load('freerunE.npz')
fr =E['xx']
from dapper import Chronology
integration_time = 4.e5
to_time = integration_time
dt = 10.0
dto = 10.0
f0 = 1.032e-4 # Coriolis parameter at 45 degrees latitude
chrono = Chronology(dt, dtObs=dto, T=to_time, BurnIn=0)
tt = chrono.tt[::chrono.dkObs]/f0/(3600*24*365) #in year
tt = tt[1:40001]
# In[]
Rmse_da = np.ones((40000,36))
Rmse_fr = np.ones((40000,36))
for i in range(40000):
mse = np.sum((da[i,0:50,:] - xx[i+1,:]) ** 2,axis=0) / 50
Rmse_da[i,:] = np.sqrt(mse.T)
mse = np.sum((fr[i,0:50,:] - xx[i+1,:]) ** 2,axis=0) / 50
Rmse_fr[i,:] = np.sqrt(mse)
# In[]
a_nx = 2
a_ny = 2
o_nx = 2
o_ny = 4
Na = a_ny*(2*a_nx+1)
No = o_ny*o_nx
Nx = 2*Na+2*No
Ns =5
o_T = 2*Na+2
o_psi = 2*Na+No+2
Ns = 5
plt.figure(figsize=(12,9))
plt.subplot(211)
plt.plot(tt,Rmse_da[:,o_psi],label = 'analysis')
#plt.plot(tt,Rmse_fr[:,o_psi],label = 'freerun')
#plt.xticks([0,9,19,27,35],fontsize = 10);
plt.legend(fontsize=10,ncol=1)
plt.title(('$\Psi$'),fontsize=10)
#plt.ylim([0,0.08]);
#plt.savefig("data RMSE.tiff", dpi=300)
plt.subplot(212)
plt.plot(tt,Rmse_da[:,o_T],label = 'analysis')
#plt.plot(tt,Rmse_fr[:,o_T],label = 'freerun')
#plt.xticks([0,9,19,27,35],fontsize = 10);
plt.legend(fontsize=10,ncol=1)
plt.title(str('$\Theta$').title(),fontsize=10)
plt.xlim([0,10])
#plt.savefig("Compare truth_DA time series.tiff", dpi=300)
plt.show()
# In[]
[xgrd,ygrd] = np.meshgrid(tt,range(36))
levels = np.linspace(-0.05, 0.05, 9)
plt.figure(figsize=(12,9))
plt.subplot(211)
plt.title('DA',fontsize=20)
plt.contourf(xgrd,ygrd,Rmse_da.T,levels=levels.round(5),cmap=plt.cm.rainbow)
plt.axhline(y=9,linestyle='--',color='black')
plt.axhline(y=19,linestyle='--',color='black')
plt.axhline(y=27,linestyle='--',color='black')
plt.yticks([9,19,27],fontsize=18)
plt.ylabel('Variable Index',fontsize=18)
plt.ylim([0,35])
plt.xlim([0,100])
plt.xticks([])
plt.colorbar()
plt.subplot(212)
plt.title('freerun',fontsize=20)
plt.contourf(xgrd,ygrd,Rmse_fr.T,levels=levels.round(5),cmap=plt.cm.rainbow)
plt.axhline(y=9,linestyle='--',color='black')
plt.axhline(y=19,linestyle='--',color='black')
plt.axhline(y=27,linestyle='--',color='black')
plt.yticks([9,19,27],fontsize=18)
plt.ylabel('Variable Index',fontsize=18)
plt.ylim([0,35])
plt.xlim([0,100])
plt.xlabel('Model Steps',fontsize=18)
plt.xlabel('Lead Times(year)',fontsize=18)
plt.colorbar()
#plt.savefig("Pcolor Compare truth_DA.tiff", dpi=300)
plt.show()
# In[]
plt.figure(figsize=(12,9))
plt.plot(Rmse_da,label = 'analysis')
plt.plot(Rmse_fr,label = 'Free Run')
plt.xlim([0,35]);
plt.axvline(x=9,linestyle='--')
plt.axvline(x=19,linestyle='--')
plt.axvline(x=27,linestyle='--')
plt.xticks([0,9,19,27,35],fontsize = 10);
plt.legend(fontsize=10,ncol=1)
plt.title(str('RMSE').title(),fontsize=10)
plt.ylim([0,0.1]);
plt.savefig("data RMSE.tiff", dpi=300)
plt.show()