-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathploteofs.py
178 lines (127 loc) · 5.27 KB
/
ploteofs.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
from netCDF4 import Dataset
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np
from eofs.standard import Eof
import constants as con
import cccmautils as cutl
# example here: http://ajdawson.github.io/eofs/examples/nao_standard.html
printtofile=True
docorr=True # EOF as correlation (otherwise covariance)
enum=1 # which EOF to retrieve and plot
# As Mori et al: 20-90N, 0-180E DJF SAT anomalies
sim='NSIDC'
field='gz50000'; ncfield='PHI'
#field='st'; ncfield='ST'
#field='pmsl'; ncfield='PMSL'
sea='DJF'
subset=True # rather than calc eof on whole globe, do a sub region
substr='morieof' # @@ make sure to also change lims below when changing this.
#substr='nh20N'
type='nh' # plot type
latlims=[20,90]
lonlims=[0,180]
#lonlims=[0,360]
limsdict = {'latlims':latlims, 'lonlims': lonlims}
if docorr: # show EOF as correlation. otherwise, covariance
cmin=-1; cmax=1
else:
cmin=''; cmax=''
fnamec,fnamep=con.build_filepathpair(sim,field)
fldc=cnc.getNCvar(fnamec,ncfield,seas=sea)
fldp=cnc.getNCvar(fnamep,ncfield,seas=sea)
lon=cnc.getNCvar(fnamec,'lon')
lat=cnc.getNCvar(fnamec,'lat')
fldca=np.squeeze(fldc-np.mean(fldc,axis=0)) # remove time mean
fldpa=np.squeeze(fldp-np.mean(fldp,axis=0))
if subset:
msk1 = con.get_t63regionmask('other',limsdict=limsdict)
fldca,msk=cutl.mask_region(fldca,lat,lon,'other',limsdict=limsdict)
fldpa,msk=cutl.mask_region(fldpa,lat,lon,'other',limsdict=limsdict)
oshape=fldca.shape # original shape
print fldca.shape
fldca1=fldca[~fldca.mask] # @@@ why do I have to invert it. so dumb.
#fldca2=fldca[msk]
fldpa1=fldpa[~fldpa.mask]
nshape=fldca1.shape # new shape
print fldca1.shape
#print fldca2.shape
# only keep region of interest?
lato=lat
lono=lon
lat=lat[np.logical_and(lat>=latlims[0],lat<=latlims[1])]
lon=lon[np.logical_and(lon>=lonlims[0],lon<=lonlims[1])]
fldca1 = fldca1.reshape( (oshape[0], len(lat),len(lon)) ) #fldcash[0]/oshape[0])) # keep time dim
fldca=fldca1
fldpa1 = fldpa1.reshape( (oshape[0], len(lat),len(lon)) )
fldpa=fldpa1
coslat = np.cos(np.deg2rad(lat)).clip(0., 1.) # why square root of cos lat? (better for EOF for some reason.)
wgts = np.sqrt(coslat)[..., np.newaxis]
solverc=Eof(fldca,weights=wgts)
if docorr:
eof1c=solverc.eofsAsCorrelation(neofs=enum)
eof1c=eof1c[enum-1,...]
else:
eof1c=solverc.eofsAsCovariance(neofs=enum)
eof1c=eof1c[enum-1,...]
eof1c=eof1c.squeeze()
eigsc=solverc.eigenvalues()
vexpc=eigsc[enum-1]/eigsc.sum() *100 # percent variance explained
fig,axs=plt.subplots(1,4)
fig.set_size_inches(12,5)
ax=axs[0]
cplt.kemmap(eof1c,lat,lon,type=type, title=sim + ' control EOF' + str(enum),axis=ax,cmin=cmin,cmax=cmax)
ax.set_ylabel(str(np.round(vexpc)))
solverp=Eof(fldpa,weights=wgts)
if docorr:
eof1p=solverp.eofsAsCorrelation(neofs=enum)
eof1p=eof1p[enum-1,...]
else:
eof1p=solverp.eofsAsCovariance(neofs=enum)
eof1p=eof1p[enum-1,...]
eof1p=eof1p.squeeze()
eigsp=solverp.eigenvalues()
vexpp=eigsp[enum-1]/eigsp.sum() *100 # percent variance explained
ax=axs[1]
cplt.kemmap(eof1p,lat,lon,type=type, title='pert EOF',axis=ax,cmin=cmin,cmax=cmax)
ax.set_ylabel(str(np.round(vexpp)))
ax=axs[2]
cplt.kemmap(eof1p-eof1c,lat,lon,type=type,title='diff of EOF',axis=ax,cmin=cmin,cmax=cmax)
# Now because I have two 'equilibrium' simulations that I want to difference,
# do I calc the eofs first for the pert and control, and then difference?
# Or do I calc the difference, remove the time mean, and calc eof? @@@@
# Also, what domain should I include in the eof calc? this is the whole globe I think.
# try taking the eof of the difference. Note it is completely different.
# This version seems more arbitrary to me, because we are differencing two
# timeseries where each year could theoretically be in any order, they just
# happen to be in this order....
# Then again, the sign of the EOF is supposedly arbitrary, so differencing 2 EOFS
# may be meaningless. For now, calc difference timeseries as pert-mean(control). Then, remove the
# time mean of that, and calc EOF
# @@@ Also, I want to actually plot the EOFs in physical units. I forget how -- check notes
# from objective analysis
#fldd=fldp-fldc
fldd=fldp-np.mean(fldc,axis=0) # this gets rid of the arbitrary timeseries-ness of it...results are very similar but not exact
fldda = np.squeeze(fldd-np.mean(fldd,axis=0))
if subset:
fldda,msk=cutl.mask_region(fldda,lato,lono,'other',limsdict=limsdict)
fldda1=fldda[~fldda.mask]
fldda1 = fldda1.reshape( (oshape[0], len(lat),len(lon)) )
fldda=fldda1
solverd=Eof(fldda,weights=wgts)
if docorr:
eof1d=solverd.eofsAsCorrelation(neofs=enum)
eof1d=eof1d[enum-1,...]
else:
eof1d=solverd.eofsAsCovariance(neofs=enum)
eof1d=eof1d[enum-1,...]
eigsd=solverd.eigenvalues()
vexpd=eigsd[enum-1]/eigsd.sum() *100 # percent variance explained
ax=axs[3]
cplt.kemmap(eof1d,lat,lon,type=type,title='Eof of diff',axis=ax,cmin=cmin,cmax=cmax)
ax.set_ylabel(str(np.round(vexpd)))
if printtofile:
if docorr:
fig.savefig(field + '_globEOF' + str(enum) + 'ascorr_' + sim + '_' + sea + '_' + substr + '_' + type + '2.pdf')
else:
fig.savefig(field + '_globEOF' + str(enum) + 'ascovar_' + sim + '_' + sea + '_' +substr + '_' + type + '2.pdf')