Global Index (short | long) | Local contents | Local Index (short | long)
Load and configure NMC.REANAL data
| This script calls | |
|---|---|
clear
cd /home/disk/hayes2/dvimont/ccm/ccm3.6/data
cd /home/disk/tao/data/reynolds/eof
filin = ['sst.mon.mean.nc'];
filin = ['ssteof5092.nc'];
nc = netcdf(filin, 'nowrite');
ts = nc{'data'}(:,:,:);
datets = nc{'time'}(:);
latts = nc{'lat'}(:);
lonts = nc{'lon'}(:);
add_offset = nc{'data'}.add_offset(:);
scale_factor = nc{'data'}.scale_factor(:);
nc = close(nc)
ts = ts * scale_factor;
ts = ts + add_offset;
[ntim, nlat, nlon] = size(ts);
% Load landmask
filin = ['landmask.nc'];
nc = netcdf(filin, 'nowrite');
lm = nc{'data'}(:);
lmlat = nc{'lat'}(:);
lmlon = nc{'lon'}(:);
nc = close(nc);
% Get CCM coordinates
cd /home/disk/hayes2/dvimont/ccm/ccm3.6/data
filin = 'cold.nc';
nc = netcdf(filin, 'nowrite');
clat = nc{'lat'}(:);
clon = nc{'lon'}(:);
nc = close(nc);
% Interpolate everything to CCM coordinates
for i = 1:ntim
cts(i,:,:) = interp2(lonts, latts, squeeze(ts(i,:,:)), clon, clat');
end
ts = cts;
clear cts
clm = interp2(lmlon, lmlat, lm, clon, clat');
clear lm
[ntim, nlat, nlon] = size(ts);
% clm = lm;
% Reshape and get rid of land points
ts = reshape(ts, ntim, nlat*nlon);
clm = reshape(clm, 1, nlat*nlon);
ocpts = find(clm == 1);
ts = ts(:,ocpts);
% Calculate ctstar and such
load ct_gr.mat
% Get SSTA's instead of SST's
%[ts, clim] = annave(ts);
% Calculate low pass filtered SSTA's
%tslp = myrunning_ave(ts, 25);
%tslp = myrunning_ave(tslp, 37);
%tshp = ts - tslp;
% Calculate ctstar and such
%ctlim = [180 270 -6 6];
%lonct = find(clon >= ctlim(1) & clon <= ctlim(2));
%latct = find(clat >= ctlim(3) & clat <= ctlim(4));
%tsnew = NaN*ones(ntim, nlat*nlon);
%tsnew(:, ocpts) = tshp;
%tsnew = reshape(tsnew, ntim, nlat, nlon);
%ctstar = squeeze(mean2(mean2(shiftdim(tsnew(:,latct,lonct),1))));
%cd /home/disk/tao/dvimont/matlab/CCM/GR/DATA
%save ct_reynolds_eof.mat ctstar
%clear tshp tslp tsnew
%[ts, clim] = annave(ts, clim, 1);
%ts(find(ts < 0.)) = zeros(size(find(ts < 0.)));
% 2. Calculate regression coefficients
%[ts, clim] = annave(ts(37:(ntim-36),:));
%ctemp = ctstar(37:(ntim-36))-mean(ctstar(37:(ntim-36)));
%ctstar = ctemp;
% Calculate GR time series
% 1. First, get rid of any points over ice
ts(find(ts < 0.)) = zeros(size(find(ts < 0.)));
% 2. Calculate regression coefficients
[ts, clim] = annave(ts(37:(ntim-36),:));
%ntim = length(37:(ntim-36));
ctemp = ctstar;
a1 = (1/(ctemp'*ctemp))*ctemp'*ts;
% 3. Calculate EOF of residual time series
[ntim, nkeep] = size(ts);
resid = ts - ctemp*a1;
tsnew = zeros(ntim, nlat*nlon);
tsnew(:,ocpts) = resid;
tsnew = reshape(tsnew, ntim, nlat, nlon);
tsnew = cosweight(tsnew, clat);
tsnew = reshape(tsnew, ntim, nlat*nlon);
resid = tsnew(:,ocpts);
% Try this by eliminating the poles
tsnew = zeros(ntim, nlat*nlon);
tsnew(:,ocpts) = resid;
tsnew = reshape(tsnew, ntim, nlat, nlon);
latmask = find(clat <= -60 | clat >= 60);
tsnew(:, latmask, :) = zeros(ntim, length(latmask), nlon);
tsnew = reshape(tsnew, ntim, nlat*nlon);
resid = tsnew(:,ocpts);
c = (resid*resid')/nkeep;
[lam, pc, per] = eof(c);
pc = ones(ntim,1) * (1./(std(pc))) .* pc;
lds = pc(:,1:28)'*resid/ntim;
% First try some plots
get_global
define_global
num = 1;
tem = zeros(1, nlat*nlon);
tem(ocpts) = lds(num,:);
tem = reshape(tem, nlat, nlon);
sd(1);
gcont(1*tem, [-5:.1:5])
dc
title('Reynolds EOF SST: EOF1 from residual time series')
xlabel('Units: 0.1 K / std(pc1)')
sd(2);
plot(1*pc(:,num));
set(gca, 'XTick',[25:60:445],'XTickLabel',[55:5:100],'YTick',[-3:3])
grid
axis([0 445 -3 3])
xlabel(['Years'])
ylabel(['STD'])
title('PC1 from residual time series')
cd /home/disk/tao/dvimont/matlab/CCM/GR/GR_Plots
% Rotate the pc's
[prot, var] = varimax(pc(:,1:28), lam(1:28), 28, 1, 'N');
frot = ((ones(ntim, 1) * (1./std(prot))) .* prot)' * resid / ntim;
frot = frot';
% Rotate the loadings
nkp = 24;
wgt = diag(lds * lds');
rlds = ((1./std(lds'))' * ones(1, size(lds, 2))) .* lds;
rlds = ((1./sqrt(wgt)) * ones(1, size(lds, 2))) .* lds;
[frot, var] = varimax(rlds', lam(1:nkp), nkp, 1, 'N');
prot = frot' * resid' / size(frot,1);
prot = prot';
ptem = prot - (ones(size(prot, 1), 1) * mean(prot));
ptem = ptem .* (ones(size(prot, 1), 1) * (1./(std(prot))));
ftem = ptem' * resid ./ ntim;
ftem = ftem';
figure(1)
num = 1;
int = 0.1;
tem = zeros(1, nlat*nlon);
tem(ocpts) = ftem(:,num);
tem = reshape(tem, nlat, nlon);
sd(1);
gcont((1 * tem), [-5:int:5])
dc
title(['Reynolds EOF SST: REOF1 from Residual SST, mode ' num2str(num) ])
xlabel(['Contour Interval ' num2str(int) ' K / std'])
sd(2);
plot(1 * ptem(:,num));
axis([ 0 409 -4 4]);
set(gca, 'XTick',[37:60:409],'XTickLabel',[65:5:100],'YTick',[-3:3])
set(gca, 'XTick',[25:60:445],'XTickLabel',[55:5:100],'YTick',[-3:3])
axis([0 409 -3 3])
xlabel(['Years'])
ylabel(['STD'])
title('Time series from REOF1 of Residual SST; Corr(pc1, rpc1) = 0.86')
grid
cd /home/disk/tao/dvimont/matlab/CCM/GR/GR_Plots