A tutorial explaining the basics of the El Niño Southern Oscillation index (ENSO) and Empirical Orthogonal Functions (EOFs). Written and narrated by Chad A. Greene of NASA Jet Propulsion Laboratory.
Requires the Climate Data Toolbox for MATLAB available here: github.com/cha...
% Load example data:
load pacific_sst.mat
% View data dimensions:
whos
% View time range:
datestr(t([1 end]))
% View temporal resolution:
mean(diff(t))
% Grid up the 1d arrays:
[Lon,Lat] = meshgrid(lon,lat);
% Calculate ENSO index:
idx = enso(sst,t,Lat,Lon);
% Plot the ENSO time series:
anomaly(t,idx)
% Label El Nino events:
text(t([396 575 792]),idx([396 575 792]),...
'El Nino!','horiz','center','vert','bot')
ylabel 'sst anomaly \circC'
% Format the date axis:
datetick
% Plot the power spectral density:
plotpsd(idx,12,'logx','lambda')
xlabel 'periodicity (years)'
set(gca,'xtick',[1:7 33])
% Calculate the mean SST:
sst_mean = mean(sst,3);
% Map the mean SST:
imagescn(lon,lat,sst_mean);
cb = colorbar;
ylabel(cb,' mean temperature {\circ}C ')
cmocean thermal % sets the colormap
% Calculate SST trend (deg/yr):
sst_trend = 365.25*trend(sst,t,3);
% Map the SST trend:
imagescn(lon,lat,10*sst_trend)
cb = colorbar;
ylabel(cb,' temperature trend {\circ}C per decade ')
cmocean('balance','pivot')
% Plot the first mode of the raw SST:
imagescn(lon,lat,eof(sst,1))
colorbar;
cmocean('balance','pivot')
title 'eof first mode'
% Mark a location of interest:
hold on
plot(lon(12),lat(10),'ks')
hold off
% Get the time series at the location of interest:
sst1 = squeeze(sst(10,12,:));
plot(t,sst1)
datetick % formats the date axis
% Remove the seasonal cycle:
sst1_ds = deseason(sst1,t);
hold on
plot(t,sst1_ds)
hold off
% Remove seasonal cycle from the data:
sst_ds = deseason(sst,t);
% Plot the first mode of the deseasoned data:
imagescn(lon,lat,eof(sst_ds,1))
colorbar
cmocean('balance','pivot')
% Remove long-term trend:
sst_ds_dt = detrend3(sst_ds);
% Recreate Fig 2a of Messie & Chavez 2011 dx.doi.org/10.1...
sst_anom_var = var(sst_ds_dt,[],3);
imagescn(lon,lat,sst_anom_var);
colorbar
title('variance of temperature')
colormap(jet) % jet is inexcusable except when recreating old plots
caxis([0 1])
% Calculate eofs:
[eof_maps,pc,expv] = eof(sst_ds_dt,6);
whos eof_maps pc expv
clf % clears the previous figure
subplot(3,2,1)
imagescn(lon,lat,eof_maps(:,:,1))
axis off
cmocean('bal','pivot')
axis image
% Plot the time series of the first mode:
subplot(3,2,2)
plot(t,pc(1,:))
axis tight
box off
datetick
subplot(3,2,3)
imagescn(lon,lat,eof_maps(:,:,2))
axis off
cmocean('bal','pivot')
axis image
subplot(3,2,6)
plot(t,pc(3,:))
axis tight
box off
datetick
subplot(3,2,5)
imagescn(lon,lat,eof_maps(:,:,3))
axis off
cmocean('bal','pivot')
axis image
subplot(3,2,4)
plot(t,pc(2,:))
axis tight
box off
datetick
sgtitle 'The first three principal components'
% Explained variance:
expv
clf
% Map the observed SST anomaly of timestep 1:
subplot(1,2,1)
h1 = imagescn(lon,lat,sst_ds_dt(:,:,1));
title 'observed sst anomaly'
cmocean bal
caxis([-1 1]*2.5)
axis image
% Map the mode 1 anomaly of timestep 1:
subplot(1,2,2)
h2 = imagescn(lon,lat,eof_maps(:,:,1)*pc(1,1));
title 'reconstructed sst anomaly'
cmocean bal
caxis([-1 1]*2.5)
axis image
sgtitle(datestr(t(1),'yyyy-mmm-dd'))
% Reconstruct an SST anomaly map:
h2.CData = eof_maps(:,:,1)*pc(1,1) + ...
eof_maps(:,:,2)*pc(2,1) + ...
eof_maps(:,:,3)*pc(3,1) + ...
eof_maps(:,:,4)*pc(4,1) + ...
eof_maps(:,:,5)*pc(5,1) + ...
eof_maps(:,:,6)*pc(6,1);
% Reconstruct sst anomalies from first 5 modes:
sst_ds_dt_r = reof(eof_maps,pc,1:5);
% Animate the two time series:
for k = 1:120
h1.CData = sst_ds_dt(:,:,k);
h2.CData = sst_ds_dt_r(:,:,k);
pause(0.1)
sgtitle(datestr(t(k),'yyyy-mmm-dd'))
end
% Tip: Check out the gif function!
% Define the Nino 3.4 box:
latrange = [-5 5];
lonrange = [-170 -120];
% Create a Nino 3.4 mask:
mask = geomask(Lat,Lon,latrange,lonrange);
% Plot the mask outline:
hold on
contour(lon,lat,mask,[.5 .5],'k')
% Area of each grid cell:
A = cdtarea(Lat,Lon);
% Area-weighted SST anomaly time series:
sst_anom_34 = local(sst_ds_dt,mask,'weight',A);
clf
anomaly(t,sst_anom_34,'color','none')
axis tight
box off
datetick
ylabel 'sst anomaly \circC'
idx = enso(sst,t,Lat,Lon);
hold on
plot(t,idx,'k')
text(t([396 575 792]),idx([396 575 792]),...
'El Nino!','horiz','center','vert','bot')
clf
plot(idx,pc(1,:),'.')
hold on
polyplot(idx,pc(1,:))
pc1_scale = trend(idx,pc(1,:))
clf
anomaly(t,idx)
hold on
plot(t,pc1_scale*pc(1,:),'color',rgb('gold'),'linewidth',2)
axis tight
datetick
text(t([396 575 792]),idx([396 575 792]),...
'El Nino!','horiz','center','vert','bot')
3 окт 2024
