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 matrixX
% 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 faster way: time as the structural dimension (1836 months), since space is bigger (180*89 = 16020 grid cells)
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
% COEFF is the loadings or EOFs - the values in the structural dimension of X or X2
% SCORE contains the PCs or "time series" - the values in the statistical dimension of X or X2
% latent contains the set of explained variances for each EOF:
explained = latent/sum(latent(:)) .*100;
% reconstruct the data
reconX = SCORE*transpose(COEFF);
% look at recon error: it is the global spatially uniform mean warming, since temporal anomalies were used
% in this permuted way where X was handed to princomp() (with space treated as the statistical dimension)
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 arrayfor 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 faster way: time as the structural dimension (1836 months), since space is bigger (180*89 = 16020 grid cells)
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% COEFF is the loadings or EOFs - the values in the structural dimension of X or X2
% SCORE contains the PCs or "time series" - the values in the statistical dimension of X or X2
% latent contains the set of explained variances for each EOF:
explained = latent/sum(latent(:)) .*100;% reconstruct the data
reconX = SCORE*transpose(COEFF);% look at recon error: it is the global spatially uniform mean warming, since temporal anomalies were used
% 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.