In [1]:
%load_ext load_style
%load_style talk.css
k-means clustering of SST anomalies in the Pacific with scikit-learn¶
In [2]:
%matplotlib inline
from matplotlib import pyplot as plt
In [3]:
import os
import numpy as np
import pandas as pd
from mpl_toolkits.basemap import Basemap as bm
In [4]:
import xray
defines a function to plot a field (must be 2D)¶
In [5]:
def plot_field(X, lat, lon, vmin, vmax, step, cmap=plt.get_cmap('jet'), ax=False, title=False, grid=False):
if not ax:
f, ax = plt.subplots(figsize=(10, (X.shape[0] / float(X.shape[1])) * 10))
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, 60), labels=[0,0,0,1])
m.drawparallels([-40, 0, 40], labels=[1,0,0,0])
m.colorbar(im)
if title:
ax.set_title(title)
reads the PCs (obtained before, see sklearn_EOF_decomposition.ipynb)¶
In [6]:
PCs = pd.read_csv('../data/EOF_ERSST_PCs.csv', index_col=0, parse_dates=True)
In [7]:
PCs.head()
Out[7]:
import the KMeans class from scikit-learn¶
In [33]:
from sklearn.cluster import KMeans
How many clusters do we want ?¶
In [34]:
nclusters = 6
instantiates the k-means class¶
In [35]:
KMeans?
In [42]:
kmeans = KMeans(init='k-means++', n_clusters=nclusters, n_init=10)
fit to the data¶
In [43]:
X = PCs.values
In [44]:
type(X)
Out[44]:
In [45]:
kmeans.fit(X)
Out[45]:
the classes (clusters) are contained in the .labels_ attributes¶
In [46]:
kmeans.labels_
Out[46]:
Calculate composite anomalies for each of the clusters¶
In [47]:
ncfname = '../data/ersst.realtime.nc'
In [48]:
dset = xray.open_dataset(ncfname)
In [49]:
dsub = dset.sel(time=slice('1980','2014'),lat=slice(-40,40), lon=slice(120,290))
In [50]:
dsub
Out[50]:
In [51]:
lat = dsub['lat'].values
lon = dsub['lon'].values
lons, lats = np.meshgrid(lon, lat)
In [52]:
labels = pd.DataFrame(kmeans.labels_, index=dsub['time'], columns=['cluster'])
In [53]:
labels.head()
Out[53]:
In [54]:
pd.unique(labels.cluster)
Out[54]:
In [55]:
c = 0
In [56]:
clussub = labels.query('cluster == {}'.format(c))
In [57]:
clussub.index
Out[57]:
In [58]:
cluster = dsub.sel(time=clussub.index).mean('time')
In [59]:
plt.imshow(cluster['anom'][0,::-1,:])
Out[59]:
In [60]:
m = bm(projection='cyl',llcrnrlat=-40,urcrnrlat=40,\
llcrnrlon=120,urcrnrlon=290,\
lat_ts=0,resolution='c')
In [61]:
f, axes = plt.subplots(nrows=3,ncols=2, figsize=(14,10))
f.subplots_adjust(hspace=0.1, wspace=0.1)
axes = axes.flatten()
for c in xrange(nclusters):
index = labels.query('cluster == {}'.format(c))
cluster = dsub.sel(time=index.index).mean('time')
ax = axes[c]
plot_field(cluster['anom'][0,:,:], lats, lons, -2, 2, 0.1, \
ax=ax, cmap=plt.get_cmap('RdBu_r'), \
title="Cluster #{}: {} months".format(c+1, len(index)))
In [ ]: