In [1]:
%load_ext load_style
%load_style talk.css
PCA of SST anomalies in the Pacific with scikit-learn¶
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
In [3]:
%matplotlib inline
load the SST data¶
The file (74 Mb) can be downloaded at ftp://ftp.niwa.co.nz/incoming/fauchereaun/ersst.realtime.nc
In [6]:
import xray
In [7]:
dset = xray.open_dataset('../data/ersst.realtime.nc')
In [8]:
dset
Out[8]:
selects the period 1980 - 2014 and the tropical Pacific domain¶
In [9]:
dsub = dset.sel(time=slice('1980','2014'), zlev=0, lat=slice(-40,40), lon=slice(120,290))
In [10]:
lat = dsub['lat'].values
lon = dsub['lon'].values
sst = dsub['anom'].values
In [11]:
sst.shape
Out[11]:
reshape in 2D (time, space)¶
In [12]:
X = np.reshape(sst, (sst.shape[0], len(lat) * len(lon)), order='F')
In [13]:
np.any(np.isnan(X))
Out[13]:
Mask the land points¶
In [14]:
type(X)
Out[14]:
In [15]:
X = ma.masked_array(X, np.isnan(X))
In [16]:
type(X)
Out[16]:
In [17]:
land = X.sum(0).mask
In [18]:
ocean = -land
keep only oceanic grid-points¶
In [19]:
X = X[:,ocean]
Standardize SST using the fit and transform methods of the sklearn.preprocessing.scaler.StandardScaler¶
In [20]:
from sklearn import preprocessing
scaler = preprocessing.StandardScaler()
In [21]:
scaler_sst = scaler.fit(X)
Once the scaler object has been 'trained' on the data, we can save it as a pickle object¶
In [22]:
from sklearn.externals import joblib
In [23]:
joblib.dump(scaler_sst, '../data/scaler_sst.pkl', compress=9)
Out[23]:
In [24]:
scaler_sst = joblib.load('../data/scaler_sst.pkl')
scales: use the transform method of the scaler object¶
In [25]:
X = scaler_sst.transform(X)
verify that mean = 0 and std = 1¶
In [26]:
X.mean()
Out[26]:
In [27]:
X.std()
Out[27]:
In [28]:
X.shape
Out[28]:
EOF decomposition¶
In [29]:
from sklearn.decomposition import pca
instantiates the PCA object¶
In [30]:
skpca = pca.PCA()
fit¶
In [31]:
skpca.fit(X)
Out[31]:
Now saves the (fitted) PCA object for reuse in operations¶
In [34]:
joblib.dump(skpca, '../data/EOF.pkl', compress=9)
Out[34]:
In [35]:
from matplotlib import style
style.use('fivethirtyeight')
In [36]:
style.available
Out[36]:
In [37]:
f, ax = plt.subplots(figsize=(6,6))
ax.plot(skpca.explained_variance_ratio_[0:10]*100)
ax.plot(skpca.explained_variance_ratio_[0:10]*100,'ro')
Out[37]:
keep number of PC sufficient to explain 70 % of the original variance¶
In [38]:
ipc = np.where(skpca.explained_variance_ratio_.cumsum() >= 0.70)[0][0]
In [39]:
ipc
Out[39]:
The Principal Components (PCs) are obtained by using the transform method of the pca object (skpca)¶
In [40]:
PCs = skpca.transform(X)
In [41]:
PCs = PCs[:,:ipc]
The Empirical Orthogonal Functions (EOFs) are contained in the components_ attribute of the pca object (skpca)¶
In [42]:
EOFs = skpca.components_
In [43]:
EOFs = EOFs[:ipc,:]
In [44]:
EOFs.shape
Out[44]:
we can the reconstruct the 2D fields (maps)¶
In [45]:
EOF_recons = np.ones((ipc, len(lat) * len(lon))) * -999.
In [46]:
for i in xrange(ipc):
EOF_recons[i,ocean] = EOFs[i,:]
In [47]:
EOF_recons = ma.masked_values(np.reshape(EOF_recons, (ipc, len(lat), len(lon)), order='F'), -999.)
In [48]:
EOF_recons.shape
Out[48]:
In [49]:
plt.imshow(EOF_recons[0,:,:], origin='lower', interpolation='nearest', aspect='auto')
plt.colorbar();
scale the Principal Components¶
In [50]:
from sklearn.preprocessing import StandardScaler
In [51]:
scaler_PCs = StandardScaler()
In [52]:
scaler_PCs.fit(PCs)
Out[52]:
In [53]:
PCs_std = scaler_PCs.transform(PCs)
In [56]:
joblib.dump(scaler_PCs, '../data/scaler_PCs.pkl')
Out[56]:
In [57]:
PCdf = pd.DataFrame(PCs_std, index = dsub['time'], \
columns = ["EOF%s" % (x) for x in xrange(1, PCs_std.shape[1] +1)])
In [58]:
PCdf.head()
Out[58]:
In [59]:
PCdf.to_csv('../data/EOF_ERSST_PCs.csv')
In [60]:
from scipy.signal import detrend
In [61]:
f, ax = plt.subplots(figsize=(12,6))
PCdf.ix[:,0].plot(ax=ax, color='k', label='PC1')
#ax.set_xlabel('period', fontsize=18)
ax.plot(PCdf.index, detrend(PCdf.ix[:,0].values), 'r', label='PC1 (trend removed)')
ax.grid('off')
ax.legend(loc=1);
In [62]:
!ipython nbconvert sklearn_EOF_decomposition.ipynb --to html