Matlab+EOF+demo+on+ERSST

This matlab code demonstrates the princomp function.
Huge .avi movie of results too big for wiki page. Sigh. Show in class. [|ncload.m] [|ERSST_EOFmovie.m] =

Here is its structure and key lines:=

% SST data grab
sst = double(ncload('/Users/bem/SST/sst.mnmean.unpacked.nc','sst')); etc. (flip, deal with bad values, grab latitude for cos(lat) weighting)

% clim12 is the annual climatology: reshape time to mm, yyyy and average over yyyy
sst12153 = reshape(sst, 89, 180, 12, 153); clim12 = nanmean(sst12153,4);

% subtract clim12 from each year to get anomalies
ssta = sst12 .*0; % right size and shape array for iy = 1:153 ssta(:,:,:,iy) = sst12(:,:,:,iy) - clim12; end

% rebin mm, yyyy back to 1 time axis for princomp call
ssta = reshape(ssta,89,180,1836);

% rebin lat,lon to 1 spatial dimension for princomp call to get a **2D data matrix** X
% want the variance budget weighted by area or cos(lat), so ssta is weighted by sqrt(cos(lat))

X = reshape(ssta.*sqrt(coslat3d),[89*180 1836]);

% The usual way: time as the 'sample' dimension, space is 'structural', need to permute, the way princomp is set up:
X2 = permute(X,[2 1]);

The Magic call to princomp:
== [COEFF,SCORE,latent,tsquare] = princomp(X); % or use X2, the "right" way. Why right? Because we have subtracted a TEMPORAL mean to get these anomalies, so we should correlate anomalies across time. If we do spatial correlation of temporal anomalies, there is a part of the data matrix X that the EOF analysis can never capture: the time series of the spatial mean. This will be clearer by example.== %%% took about a minute on my desktop machine

% latent contains the set of explained variances for each EOF:
explained = latent/sum(latent) .*100;

% reconstruct the data
reconX = SCORE*transpose(COEFF);

% in this permuted way where X was handed to princomp (with //space// treated as the statistical dimension)
err = X-reconX; err3d = reshape(err,89,180,1836); plot(mean(err,1))

% recon to truncation 5 by zeroing out COEFFs
COEFF2 = COEFF; COEFF2(:,5+1:end) = 0; recon5 = SCORE*transpose(COEFF2);

Structure of individual modes:

Huge .avi movie of results too big for wiki page. Sigh. Show in class.