(1) Recreate Figure 1 from Horel and Wallace (1981) using the NCEP reanalysis data.
In [2]:
import pandas as pd
import numpy as np
from numpy import ma
from matplotlib import pyplot as plt
from mpl_toolkits.basemap import Basemap as bm
import xray
%matplotlib inline
In [3]:
sstdata=xray.open_dataset('sst_seas_ncep.nc')
sstdata
Out[3]:
In [4]:
#find time series at point nearest to 10N and 95W
tser=sstdata.sel(lat=[10.],lon=[360.-95.],method='nearest')
tser
plt.figure()
plt.plot(tser['time'].values,tser['skt'].values.ravel())
#sst_resample
Out[4]:
In [5]:
ls=xray.open_dataset('lsmask.19294.nc',decode_times=False)
ls['lsmask']
Out[5]:
In [6]:
plt.figure()
plt.imshow(ls['lsmask'].values[0,:,:])
plt.colorbar
Out[6]:
In [7]:
mask3d=ls['lsmask'].values.repeat(272,0)
np.shape(mask3d)
Out[7]:
In [8]:
sstdata
Out[8]:
In [9]:
sstdata['skt'].values=np.ma.masked_where(mask3d == -1, sstdata['skt'].values)
plt.figure()
plt.pcolormesh(sstdata['lon'].values,sstdata['lat'].values,sstdata['skt'].values[0,:,:])
plt.clim([-20,40])
plt.colorbar()
Out[9]:
In [10]:
#calculate seasonal climatology
sstseas=sstdata.groupby('time.season').mean(dim='time')
for i in np.arange(4):
plt.figure()
plt.pcolormesh(sstseas['lon'].values,sstseas['lat'].values,sstseas['skt'].values[i,:,:])
plt.clim([-20,40])
plt.title(sstseas['season'].values[i])
plt.colorbar()
In [11]:
import numpy as np
In [16]:
fixed_sstseas=xray.concat([sstseas.sel(season='DJF'),sstseas.sel(season='MAM'),sstseas.sel(season='JJA'),sstseas.sel(season='SON')],'season')
#resize (tile ) to be same shape as data
meansstseas=np.tile(fixed_sstseas['skt'].values,(68,1,1))
anomsst=sstdata['skt'].values-meansstseas
plt.pcolormesh(sstdata['lon'].values,sstdata['lat'].values,anomsst[4,:,:])
plt.clim([-10,10])
plt.colorbar()
Out[16]:
In [18]:
def seas_departure(inputarr,field):
seas=inputarr.groupby('time.season').mean(dim='time')
fixed_seas=xray.concat([seas.sel(season='DJF'),seas.sel(season='MAM'),\
seas.sel(season='JJA'),seas.sel(season='SON')], 'season')
#calculates the seasonal climatology but np.tile gives you back a numpy array
mean_seas=np.tile(fixed_seas[field].values,(68,1,1))
#put it back into an xray dataset
#foo = xray.DataArray(data, coords=[times, locs], dims=['time', 'space'])
mean_seas=xray.DataArray(mean_seas,coords=[inputarr['time'].values,inputarr['lat'].values,\
inputarr['lon'].values],dims=['time','lat','lon'])
anom = inputarr[field].values-mean_seas
return anom
anomsst=seas_departure(sstdata,'skt')
anomsst
Out[18]:
In [19]:
dsub = sstdata.sel(lat=slice(40,-40), lon=slice(120,290))
lat = dsub['lat'].values
lon = dsub['lon'].values
#calculate sst climatology
meansst=dsub['skt'].mean(axis=0)
plt.figure()
plt.pcolormesh(meansst.values)
plt.clim([-20,40])
plt.colorbar
#resize to be the same shape as the data
meansst=meansst.values.reshape(1,np.shape(meansst)[0],np.shape(meansst)[1]).repeat(272,0)
#anomalies
anomsst=dsub['skt'].values-meansst
lons,lats = np.meshgrid(lon,lat)
plt.figure()
plt.pcolormesh(anomsst[2,:,:]) #summer in NH
plt.clim([-5,5])
plt.colorbar
Out[19]:
In [20]:
X = np.reshape(anomsst, (anomsst.shape[0], len(lat) * len(lon)), order='F')
X.shape
Out[20]:
In [21]:
X = ma.masked_array(X, np.isnan(X))
type(X)
land = X.sum(0).mask
ocean=-land
X = X[:,ocean]
X.shape
Out[21]:
In [22]:
from sklearn import preprocessing
scaler = preprocessing.StandardScaler()
scaler_sst = scaler.fit(X)
X = scaler_sst.transform(X)
from sklearn.decomposition import pca
skpca = pca.PCA()
skpca.fit(X)
Out[22]:
In [23]:
f, ax = plt.subplots(figsize=(5,5))
ax.plot(skpca.explained_variance_ratio_[0:10]*100)
ax.plot(skpca.explained_variance_ratio_[0:10]*100,'ro')
ax.set_title("% of variance explained", fontsize=14)
ax.grid()
In [24]:
PCs = skpca.transform(X)
EOFs = skpca.components_
ipc=272
EOF_recons = np.ones((ipc, len(lat) * len(lon))) * -999.
for i in range(ipc):
EOF_recons[i,ocean] = EOFs[i,:]
EOF_recons = ma.masked_values(np.reshape(EOF_recons, (ipc, len(lat), len(lon)), order='F'), -999.)
EOF_recons.shape
Out[24]:
In [25]:
def plot_field(m, X, lats, lons, vmin, vmax, step, cmap=plt.get_cmap('jet'), \
ax=False, title=False, grid=False):
if not ax:
f, ax = plt.subplots(figsize=(8, (X.shape[0] / float(X.shape[1])) * 8))
m.ax = ax
im = m.contourf(lons, lats, X, np.arange(vmin, vmax+step, step), \
latlon=True, cmap=cmap, extend='both', ax=ax)
m.drawcoastlines()
if grid:
m.drawmeridians(np.arange(0, 360, 30), labels=[0,0,0,1])
m.drawparallels(np.arange(-80, 80, 20), labels=[1,0,0,0])
m.colorbar(im)
if title:
ax.set_title(title)
m = bm(projection='cyl',llcrnrlat=lat.min(),urcrnrlat=lat.max(),\
llcrnrlon=lon.min(),urcrnrlon=lon.max(),\
lat_ts=0,resolution='c')
plot_field(m, EOF_recons[1,:,:], lats,lons, -.1 ,.1, .01, grid =True, cmap=plt.get_cmap('RdBu_r'))
In [26]:
from sklearn.preprocessing import StandardScaler
scaler_PCs = StandardScaler()
scaler_PCs.fit(PCs)
PCs_std = scaler_PCs.transform(PCs)
PCdf = pd.DataFrame(PCs_std, index = dsub['time'], \
columns = ["EOF%s" % (x) for x in range(1, PCs_std.shape[1] +1)])
from scipy.signal import detrend
f, ax = plt.subplots(figsize=(10,4))
PCdf.ix[:,1].plot(ax=ax, color='k', label='PC1')
ax.axhline(0, c='0.8')
#ax.set_xlabel('period', fontsize=18)
ax.plot(PCdf.index, detrend(PCdf.ix[:,1].values), 'r', label='PC1 (trend removed)')
ax.grid('off')
ax.legend(loc=1);
(2) Create normalized index values (mean=0, sd=1) of Sea Level Pressure between Tahiti and Darwin similar to Fig. 4. Invert the index (multiply by -1, as in the paper). Plot the time series.
In [28]:
slp=xray.open_dataset('slp_seas_ncep.nc')
slp
tahitimslp = slp.sel(lat=-15., lon=360-149.,method='nearest')
tahitimslp
darwinmslp = slp.sel(lat=-12., lon=131.,method='nearest')
darwinmslp
plt.figure()
plt.plot(tahitimslp['time'].values,-1*(tahitimslp['slp'].values-darwinmslp['slp'].values))
plt.title('-1*Tahiti-Darwin mslp')
plt.ylabel('MSLP (hPa)')
plt.xlabel('Time (years)')
#calculate normalized index, with mean=0 and std=1
diff=-1*(tahitimslp['slp']-darwinmslp['slp'])
diff=(diff-np.mean(diff))
diff=diff/np.std(diff)
diff
plt.figure()
plt.plot(diff['time'].values,diff)
plt.title('normalized Tahiti-Darwin mslp index')
plt.ylabel('index')
plt.xlabel('Time (years)')
Out[28]:
(3) Create normalized index values of 200 hPa height at the stations listed in Fig 6. Plot the time series.
In [29]:
h200=xray.open_dataset('h200_seas_ncep.nc')
#west(lon) and south(lat) negative
sites= ['NAI' ,'DAR' , 'NAN' ,'HIL' , 'SAJ']
lat=np.array([-1.2, -12.4, -17.5, 19.7, 18.0])
lon=np.array([36.8, 130.0, 177.0, -155.0, -66.0])
def grab_h200(h200,lat,lon):
return h200.sel(lat=lat,lon=360-lon, method='nearest')
all_h200={}
for i in np.arange(len(lat)):
all_h200[sites[i]]=grab_h200(h200,lat[i],lon[i])
z200index=0.2*all_h200['NAI']+0.23*all_h200['DAR']+0.23*all_h200['NAN']+0.25*all_h200['HIL']+0.27*all_h200['SAJ']
plt.plot(z200index['time'].values,z200index['hgt'].values)
def standardize(var):
var=var-np.mean(var)
return var/np.std(var)
z200index_std=standardize(z200index)
(4) Compute the PNA and WP indices as on p. 823 of the paper. Plot the time series.
In [ ]:
(5) Compute standardized rainfall values at the locations in Fig. 3. Plot the time series.
In [ ]:
(6) Reproduce Fig 8 with the top 4 and bottom 4 values of PNA index (that you find in part 4), except contour the values you have rather than plot points as in the paper.
In [30]:
#200 hPa - 700 hPa correlation map
z200index_std['hgt']
h700=xray.open_dataset('h700_seas_ncep.nc')
print(np.shape(h700['hgt'].values))
def spatial_corr(index,mapvals):
sz=np.shape(mapvals)
outmap=np.zeros([sz[1],sz[2]])
for i in np.arange(sz[1]):
for j in np.arange(sz[2]):
maptser=np.squeeze(mapvals[:,i,j])
#print(np.shape(maptser))
#print(np.shape(index))
df_bis = pd.DataFrame({'maptser': maptser,
'index': np.squeeze(index)})
outmap[i,j]=df_bis.corr().values[1,0]
return outmap
In [31]:
outmap=spatial_corr(z200index_std['hgt'].values,h700['hgt'].values)
In [32]:
from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
# setup north polar stereographic basemap.
# The longitude lon_0 is at 6-o'clock, and the
# latitude circle boundinglat is tangent to the edge
# of the map at lon_0. Default value of lat_ts
# (latitude of true scale) is pole.
m = Basemap(lon_0=180)
m.pcolormesh(h700['lon'].values,h700['lat'].values,outmap)
m.drawcoastlines()
plt.colorbar()
Out[32]:
(7) Recreate Figure 9 using your data.
In [ ]: