Matlab Implementation Example of Empirical Orthogonal Decomposition EOF

In geosciences, PCA and EOF are usually used for signal extraction. Separating the spatiotemporal variation characteristics of geographic elements from complex spatiotemporal data is a prerequisite for geoscience signal analysis. In essence, there is no difference between PCA and EOF, except that: EOF is a spatial feature vector, also known as a spatial mode, which reflects the spatial distribution characteristics of the element field to a certain extent; PC (principal component) corresponds to time changes, also known as time Coefficient, reflecting the change in weight of the corresponding spatial mode (EOF) over time. In short, they are two important elements in the decomposition process using empirical orthogonality.

Below I skip the introduction of the basic concepts, and directly give a sample code for Matlab to perform EOF decomposition:

1. Load data

% Load sample data:
load pacific_sst.mat
% Calculate the first EOF of sea surface temperatures and its
% principal component time series:
[eofmap,pc] = eof(sst,1);
% Plot the first EOF map:
imagescn(lon,lat,eofmap);
axis xy image off
% Optional: Use a cmocean colormap:
cmocean('delta','pivot',0)

The result is the first EOF of the sea temperature dataset

But since we didn’t remove the seasonal period, the first EOF mainly represents the seasonal change. As evidence that the above pattern is related to the seasonal cycle, look at the corresponding principal component time series.

figure
anomaly(t,pc)
axis tight
xlim([datenum('jan 1, 1990') datenum('jan 1, 1995')])
datetick('x','keeplimits')

2. Calculate the average ocean temperature data and display the average temperature for many years

figure
imagescn(lon,lat,mean(sst,3));
axis xy off
cb = colorbar;
ylabel(cb,' mean temperature {\circ}C ')
cmocean thermal

3. Calculate the linear change trend from 1950 to 2016, making sure to multiply the trend by 10*365.25 to convert degrees per day to degrees per decade:

In order to remove the influence of seasonal items and trend items on abnormal signals, we next detrend and decycle the data

sst = detrend3(sst,t);
sst = deseason(sst,t);

So now our sea temperature data has been detrended, the mean has been removed, and the seasonal cycle has been removed. All that’s left in SST are anomalies — things that vary, but not long-term trends or short-term annual cycles. Here is the residual variance for our SST anomaly dataset:

4. Then calculate the first three modes of EOF

figure
subsubplot(3,1,1)
plot(t,pc(1,:))
box off
axis tight
ylabel 'pc1'
title 'The first three principal components'

subsubplot(3,1,2)
plot(t,pc(2,:))
box off
axis tight
set(gca,'yaxislocation','right')
ylabel 'pc2'

subsubplot(3,1,3)
plot(t,pc(3,:))
box off
axis tight
ylabel 'pc3'
datetick('x','keeplimits')

5. El Ni?o Southern Oscillation (ENSO) time series, sure enough, the strongest El Ni?o events on record occurred in 1982-1983, 1997-1998 and 2014-2016.

figure('pos',[100 100 600 250])
anomaly(t,-pc(1,:),'topcolor',rgb('bubblegum'),...
   'bottomcolor',rgb('periwinkle blue')) % First principal component is enso
box off
axis tight
datetick('x','keeplimits')
text([724316 729713 736290],[.95 .99 .81],'El Nino','horiz','center')

Sometimes we hear that the characteristic frequency of El Ni?o is every 5 years, or 5 to 7 years, sometimes you hear it is every 2 to 7 years. This is hard to see in a time series, so we plot the first principal component in the frequency domain with plotpsd, specifying a sampling frequency of 12 samples per year, plotted on the log x-axis, with x-values in lambda(years), instead of frequency:

figure
plotpsd(pc(1,:),12,'logx','lambda')
xlabel 'periodicity (years)'
set(gca,'xtick',[1:7 33])

Below is the spatial distribution of the first 6 modalities:

s = [-1 1 -1 1 -1 1]; % (sign multiplier to match Messie and Chavez 2011)

figure('pos',[100 100 500 700])
for k = 1:6
   subplot(3,2,k)
   imagescn(lon,lat,eof_maps(:,:,k)*s(k));
   axis xy off
   title(['Mode ',num2str(k),' (',num2str(expv(k),'%0.1f'),'%)'])
   caxis([-2 2])
end

colormap jet

At any given time, a snapshot of the ENSO-related SST anomaly can be obtained by plotting Mode 1 as shown above, multiplied by its corresponding principal component (vector pc(1,:))). Likewise, you can get a picture of global sea surface temperature anomalies for a given time by adding up all EOF maps and multiplying by the corresponding principal components at that time. In this way, we can build an increasingly complete movie of SST anomalies as we include more and more modes of variability.

For example, a map of SST anomalies associated with the first three modes of variability at a specific time can be obtained by summing the EOF maps for each mode, multiplying by the corresponding pc values. Calculate the sum of the 90s like this:

% Indices of start and end dates for the movie:
startind = find(t>=datenum('jan 1, 1990'),1,'first');
endind = find(t<=datenum('dec 31, 1999'),1,'last');
% A map of SST anomalies from first three modes at start:
map = eof_maps(:,:,1)*pc(1,startind) + ... % Mode 1, Jan 1990
      eof_maps(:,:,2)*pc(2,startind) + ... % Mode 2, Jan 1990
      eof_maps(:,:,3)*pc(3,startind); % Mode 3, Jan 1990
sst_f = reof(eof_maps,pc,1:3);
ind_1990s = 481:3:600; % (every third value to cut down on gif size)
figure
h = imagescn(sst_f(:,:,ind_1990s(1)));
caxis([-2 2])
cmocean balance
cb = colorbar;
ylabel(cb,'temperature anomaly (modes 1-3)')
title(datestr(t(ind_1990s(1)),'yyyy'))
gif('SSTs_1990s.gif','frame',gcf) % writes the first frame
for k = 2:length(ind_1990s)
   h.CData = sst_f(:,:,ind_1990s(k));
   title(datestr(t(ind_1990s(k)),'yyyy'))
   gif % adds this frame to the gif
end

References:

  1. The principle of principal component analysis PCA and empirical orthogonal function analysis EOF (easy-to-understand explanation) – a99h’s blog – CSDN blog

  2. Messié, Monique, and Francisco Chavez. “Global modes of sea surface temperature variability in relation to regional climate indices.” Journal of Climate 24.16 (2011): 4314-4331. doi:10.1175/2011JCLI3941.1.

  3. Rayner, N. A., Parker, D. E., Horton, E. B., Folland, C. K., Alexander, L. V., Rowell, D. P., Kent, E. C., Kaplan, A. (2003). Global analyzes of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century. J. Geophys. Res. Vol. 108, No. D14, 4407 doi:10.1029/2002JD002670.

  4. Thyng, K.M., C.A. Greene, R.D. Hetland, H.M. Zimmerle, and S.F. DiMarco. 2016. True colors of oceanography: Guidelines for effective and accurate colormap selection. Oceanography 29(3):9-13, doi:10.5670/oceanog.2016.

Matlab program and data download:

Link: https://pan.baidu.com/s/1IrlhONl8V9ubuykuy7yq_g?pwd=ucas

Extraction code: ucas

Welcome to like and collect!