In [1]:
%load_ext load_style
%load_style talk.css
In statistics and signal processing, the method of empirical orthogonal function (EOF) analysis is a decomposition of a signal or data set in terms of orthogonal basis functions which are determined from the data. It is similar to performing a principal components analysis on the data, except that the EOF method finds both temporal projections and spatial patterns. The term is also interchangeable with the geographically weighted PCAs in geophysics (https://en.wikipedia.org/wiki/Empirical_orthogonal_functions).
The spatial patterns are the EOFs, and can be thought of as basis functions in terms of variance. The associated temporal projections are the pricipal components (PCs) and are the temporal coefficients of the EOF patterns.
This notebook is inspired by the blog presented by yiboj, where you can download the sample code and sample data. For convinience, yearly data have been converted into a single file using CDO. However, it should be noted this notebook simplifies the original source code through some advanced syntaxes of NumPy and fixs a bug of filling zeros over land.
The SST data is from from AVHRR Level 4 dataset in central pacific ocean from 1982 to 2000.
In [2]:
% matplotlib inline
import numpy as np
from scipy import signal
import numpy.polynomial.polynomial as poly
from netCDF4 import Dataset
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from eofs.standard import Eof
In [3]:
startY = 1982
endY = 2000
ndy = 36
infile = 'data\eof_data\sst.sw.AVHRR.l4.1982.2000.nc'
ncin = Dataset(infile, 'r')
sst_raw = ncin.variables['analysed_sst'][:]
lons = ncin.variables['lon'][:]
lats = ncin.variables['lat'][:]
ncin.close()
nt,nlat,nlon = sst_raw.shape
ny = nt/ndy
mask = sst_raw.mask
In [4]:
sst_detrend = np.empty(sst_raw.shape)
sst_coeffs = np.empty((2, nlat, nlon))
sst_detrend[:,:,:] = np.nan
In [5]:
x = np.linspace(1,nt,nt)
for i in range(0, nlat):
for j in range(0,nlon):
ytemp = np.copy(sst_raw[:,i,j])
y = sst_raw[:,i,j]
b = ~np.isnan(y)
coefs = poly.polyfit(x[b], y[b], 1)
sst_coeffs[0,i,j] = coefs[0]
sst_coeffs[1,i,j] = coefs[1]
ffit = poly.polyval(x[b], coefs)
sst_detrend[b,i,j] = y[b] - ffit
In [6]:
sst_all = sst_detrend.reshape((ndy,ny,nlat,nlon), order='F').transpose((1,0,2,3)) # year, 36, lat, lon
In [7]:
sst_season = np.mean(sst_all, axis=0)
In [8]:
sst_diff = sst_all - sst_season
sst_diff = np.ma.masked_array(sst_diff, mask=mask) # have to do this, or fill in zeros in sst_diff.
In [9]:
sst_final = sst_diff.transpose((1,0,2,3)).reshape((ndy*ny,nlat,nlon), order='F')
In [10]:
coslat = np.cos(np.deg2rad(lats))
wgts = np.sqrt(coslat)[..., np.newaxis]
solver = Eof(sst_final, weights=wgts)
print(coslat.shape)
print(wgts.shape)
In [11]:
eof1 = solver.eofs(neofs=10)
pc1 = solver.pcs(npcs=10, pcscaling=0)
varfrac = solver.varianceFraction()
lambdas = solver.eigenvalues()
In [12]:
parallels = np.arange(-90,90,10.)
meridians = np.arange(-180,180,20)
for i in range(0, 3):
fig = plt.figure(figsize=(9,7))
plt.subplot(211)
#ax=fig.add_axes([0.1,0.1,0.8,0.8])
m = Basemap(projection='cyl', llcrnrlon=min(lons), llcrnrlat=min(lats),
urcrnrlon=max(lons), urcrnrlat=max(lats))
x, y = m(*np.meshgrid(lons, lats))
clevs = np.linspace(np.min(eof1[i,:,:].squeeze()), np.max(eof1[i,:,:].squeeze()), 11)
cs = m.contourf(x, y, eof1[i,:,:].squeeze(), clevs, cmap=plt.cm.RdBu_r)
m.drawcoastlines()
#m.fillcontinents(color='#000000',lake_color='#99ffff')
m.drawparallels(parallels,labels=[1,0,0,0])
m.drawmeridians(meridians,labels=[1,0,0,1])
#cb = plt.colorbar(cs, orientation='horizontal')
cb = m.colorbar(cs, 'right', size='5%', pad='2%')
cb.set_label('EOF', fontsize=12)
plt.title('EOF ' + str(i+1), fontsize=16)
plt.subplot(212)
days = [startY+(x*10+1)/365.0 for x in range(0, nt)]
plt.plot(days, pc1[:,i], linewidth=2)
plt.xticks(range(startY, endY), rotation='vertical')
plt.axhline(0, color='k')
plt.xlabel('Year')
plt.ylabel('PC Amplitude')
plt.xlim(startY, endY)
plt.ylim(np.min(pc1.squeeze()), np.max(pc1.squeeze()))
plt.xticks(range(startY, endY))
plt.tight_layout()
In [13]:
plt.figure(figsize=(11,6))
eof_num = range(1, 16)
plt.plot(eof_num, varfrac[0:15], linewidth=2)
plt.plot(eof_num, varfrac[0:15], linestyle='None', marker="o", color='r', markersize=8)
plt.axhline(0, color='k')
plt.xticks(range(1, 16))
plt.title('Fraction of the total variance represented by each EOF')
plt.xlabel('EOF #')
plt.ylabel('Variance Fraction')
plt.xlim(1, 15)
plt.ylim(np.min(varfrac), np.max(varfrac)+0.01)
Out[13]:
The spatial pattern of this EOF is the canonical El Nino pattern, and the associated time series shows large peaks and troughs for well-known El Nino and La Nina events.
http://unidata.github.io/netcdf4-python/
John D. Hunter. Matplotlib: A 2D Graphics Environment, Computing in Science & Engineering, 9, 90-95 (2007), DOI:10.1109/MCSE.2007.55
Stéfan van der Walt, S. Chris Colbert and Gaël Varoquaux. The NumPy Array: A Structure for Efficient Numerical Computation, Computing in Science & Engineering, 13, 22-30 (2011), DOI:10.1109/MCSE.2011.37
In [ ]: