TimeSeries_smblike.py 4.82 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
# Plot time series plot of Figure 13 in the paper 

import numpy as np
import numpy.ma as ma
#import Nio,Ngl,os
import matplotlib.pyplot as plt
plt.switch_backend('agg')

from netCDF4 import Dataset

# ---------------- Data source in rhea, the two files save the same data ------------------------
12
13
14
15
input_file11 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/tseries/20180215.DECKv1b_H1.ne30_oEC.edison.clm2.aavg.h0.yrly.RAIN.nc'
input_file12 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/tseries/20180215.DECKv1b_H1.ne30_oEC.edison.clm2.aavg.h0.yrly.SNOW.nc'
input_file13 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/tseries/20180215.DECKv1b_H1.ne30_oEC.edison.clm2.aavg.h0.yrly.QRUNOFF.nc'
input_file14 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/tseries/20180215.DECKv1b_H1.ne30_oEC.edison.clm2.aavg.h0.yrly.QSOIL.nc'
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
input_file2 = '/lustre/atlas1/cli115/world-shared/4ue/racmo23_GRN_monthly/tseries/racmo23_GRN_monthly.smb.ann_tseries_aavg.nc'
# --------------------------------------------------------------

# read input_file11 
# time from X-X, both time and acab have 156 data
ncid11 = Dataset(input_file11, mode='r')
time = ncid11.variables['mcdate'][:]
rain = ncid11.variables['RAIN'][:] #m/yr
ncid11.set_auto_scale(True) # put the scale-factor in the time-series plot

#time = (time-time[1])/365 + 1851 # time from 1851-2006
timeB = [x for x in range(1, 1+len(time))] 
timeB = [x + 1960 for x in timeB]
#print "timeB:", timeB

# read input_file12 
# time from X-X, both time and acab have 156 data
ncid12 = Dataset(input_file12, mode='r')
time = ncid12.variables['time'][:]
snow = ncid12.variables['SNOW'][:] #m/yr
ncid12.set_auto_scale(True) # put the scale-factor in the time-series plot

# read input_file13 
# time from X-X, both time and acab have 156 data
ncid13 = Dataset(input_file13, mode='r')
qrun = ncid13.variables['QRUNOFF'][:] #m/yr
ncid13.set_auto_scale(True) # put the scale-factor in the time-series plot

# read input_file14 
# time from X-X, both time and acab have 156 data
ncid14 = Dataset(input_file14, mode='r')
qsoil = ncid14.variables['QSOIL'][:] #m/yr
ncid14.set_auto_scale(True) # put the scale-factor in the time-series plot

# read input_file2 
ncid2 = Dataset(input_file2, mode='r')
time2 = ncid2.variables['time'][:]
smb = ncid2.variables['smb'][:]  # mmWE
#ncid2.set_auto_scale(True) # put the scale-factor in the time-series plot

#smb = smb #unit transfer
time2 = (time2+6.5)/12 + 1957 # add months to make the full year, divide my months to get years, add to get 1958
#print "time2:", time2

qice = (rain + snow - qrun - qsoil)*3.1567e7 # convert to kg/m^2/yr

maxqice=(np.max(qice))
minqice=(np.min(qice))
maxtime=(np.max(timeB))
mintime=(np.min(timeB))
maxsmb=(np.max(smb))
minsmb=(np.min(smb))
maxtime2=(np.max(time2))
mintime2=(np.min(time2))

print "Min QICE E3SM:", minqice
print "Max QICE E3SM:", maxqice
#print "Max time:", maxtime
#print "Min time:", mintime
print "Min smb RACMO:", minsmb
print "Max smb RACMO:", maxsmb
#print "Max time2:", maxtime2
#print "Min time2:", mintime2

##------- Time Series PLOT --------
fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(111)

ax.plot(timeB,qice,'r-')
ax.plot(time2,smb,'b-')

ax.set_ylim([-100, 400])
ax.set_xlim([1955, 2015])
ax.xaxis.set_ticks(np.arange(1950,2010,10))
#ax.set_xticklabels(['1950','1960','1970','1980','1990','2000','2010'])
ax.set_ylabel('SMB (kg/m2/yr)',fontsize=14)
ax.set_xlabel('(years)',fontsize=14)
##ax.text(-12,5,'T-2m, JJA',fontsize=18)
##ax.legend([s1,s2],['>20%','>99%'],loc='lower right')

96
plt.savefig('/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/plots/TimeSeries_smblike.ps')
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
#plt.savefig('TimeSeries_SMB.ps')
plt.close()

##------ Box plot of statistics ---------
## If I have 3 sets of data, put them a matrix with 3 columns,
## then boxplot(matrix) in one figure data=[a1,a2,a3]

fig = plt.figure(figsize=(4.8,4))

#--- the whisker shows [Q1-1.5IQR, Q3+1.5IQR] without outliners
#ax = fig.add_subplot(121)
#ax.boxplot(qice,0,'')
#ax.set_ylim([0,0])
#ax.set_ylabel('qice(kg/m2/yr)',fontsize=14)
#ax.set_xticklabels('')

#--- the whisker shows [Q1-1.5IQR, Q3+1.5IQR] with outliners
#--- where the Q1 is 25% percentile and Q3 is 75% percential, IQR=Q3-Q1
#--- the outlier is data point located outside the whiskers.
ax = fig.add_subplot(121)
ax.boxplot(qice,0,'gd')
ax.set_ylabel('E3SM SMB(km/m2/yr)',fontsize=14)
ax.set_ylim([0,500])
ax.set_xticklabels('')
ax = fig.add_subplot(122)
ax.boxplot(smb,0,'gd')
ax.set_ylabel('RACMO SMB(km/m2/yr)',fontsize=14)
ax.set_ylim([0,500])
ax.set_xticklabels('')

127
plt.savefig('/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/plots/BoxPlot_smblike.ps')
128
129
130
131
#plt.savefig('BoxPlot_SMB.ps')
plt.close()