Global Index (short | long) | Local contents | Local Index (short | long)
cd /home/disk/tao/data/coads/coads1 filin = 'sst.mean.nc';
| This script calls | |
|---|---|
clear
cd /home/disk/tao/data/coads/sstanom4by6
filin = 'sstcoadsanom4by6.18541993.nc'
lim1 = 553;
%lim2 = 1512;
lim2 = 1680;
varnam = 'data';
nc = netcdf(filin, 'nowrite');
ts = nc{varnam}(lim1:lim2,:,:);
datets = nc{'time'}(lim1:lim2);
lat = nc{'lat'}(:);
lon = nc{'lon'}(:);
add_offset = nc{varnam}.add_offset(:);
scale_factor = nc{varnam}.scale_factor(:);
mv = nc{varnam}.missing_value(:);
nc = close(nc)
ts = ts * scale_factor;
ts = ts + add_offset;
mv = mv * scale_factor;
mv = mv + add_offset;
cd /home/disk/tao/data/coads/coads1a
filin = 'sst.mean.8095.nc';
nc = netcdf(filin, 'nowrite');
ts2 = nc{var}(:,:,:);
add_offset = nc{var}.add_offset(:);
scale_factor = nc{var}.scale_factor(:);
mv2 = nc{var}.missing_value(:);
nc = close(nc);
ts2 = ts2 * scale_factor;
ts2 = ts2 + add_offset;
mv2 = mv2 * scale_factor;
mv2 = mv2 + add_offset;
ts2(find(ts2 == mv2)) = mv * ones(size(find(ts2 == mv2)));
ts = [ts; ts2];
clear ts2
[ntim, nlat, nlon] = size(ts);
ts = reshape(ts, ntim, nlat*nlon);
tsave = mean(ts);
kp = find(tsave ~= mv);
ts = ts(:, kp);
kp2 = [];
for i = 1:length(kp);
sznan = length(find(ts(:, i) == mv));
if sznan < 1002;
kp2 = [kp2, i];
end;
end;
kp3 = find(ismember(kp, kp(kp2)));
ts = ts(:, kp3);
kp = kp(kp3);
ts(find(ts == mv)) = NaN * ones(size(find(ts == mv)));
[ts, clim] = annave(ts);
% Find CT*:
getlp = 0;
if getlp;
tslp = myrunning_ave(ts, 25);
tslp = myrunning_ave(tslp, 37);
tshp = ts - tslp;
cd /home/disk/hayes2/dvimont/data
save coads_2x2_sst.mat tshp tslp ts;
else
cd /home/disk/hayes2/dvimont/data
load coads_2x2_sst.mat
ntim = size(ts, 1);
nlat = length(lat);
nlon = length(lon);
end
%kptime = 37:(ntim - 36)
%kptime = 1:ntim;
kptime = 1:612;
ntim = 612;
clear tslp
%[tshp, clim] = annave(tshp);
tsnew = NaN * ones(length(kptime), nlat*nlon);
tsnew(:, kp) = ts(kptime, :);
tsnew = reshape(tsnew, length(kptime), nlat, nlon);
ctlim = [180 270 -6 6];
[xk, yk] = keep_var(ctlim, lon, lat);
ctstar = squeeze(mean2(mean2(shiftdim(tsnew(kptime, yk, xk), 1))));
ct = (ctstar - mean2(ctstar)) / std(ctstar(find(~isnan(ctstar))));
ctpat = regress_nan(ctstar, ts(kptime, :));
% Highpass filter:
ctstar = ct - myrunning_ave(myrunning_ave(ct, 25),37);
ctstar = (ctstar - mean2(ctstar)) / std(ctstar(find(~isnan(ctstar))));
cthpat = regress_nan(ctstar, ts(kptime, :));
tem = NaN * ones(1, nlat*nlon);
tem(kp) = cthpat;
tem = reshape(tem, nlat, nlon);
get_global
default_global
FRAME = [min(lon) (min(lon) + 360) -90 90];
figure(1)
sd(1)
gcont(tem, [-1:.2:1]);
dc
sd(2)
plot(ctstar)
set(gca, 'XTick', [1:120:1153], 'XTickLabel', [1900:10:2000])
axis([0 ntim+1 -3 4])
grid
% Now look at EOF's:
kptime = 1:500;
clear tshp
[ts, dum] = annave(ts(kptime, :));
ntim = size(ts, 1);
tsnew = NaN * ones(ntim, nlat*nlon);
tsnew(:, kp) = ts;
tsnew = reshape(tsnew, ntim, nlat, nlon);
tsnew = cosweight(tsnew, lat);
tsnew = reshape(tsnew, ntim, nlat*nlon);
tsnew = tsnew(:, kp);
%c = covar_nan(tsnew', tsnew');
% Try a different approach:
n_not_nan = zeros(ntim);
tsc = zeros(size(tsnew));
for i = 1:ntim;
kpc = find(~isnan(tsnew(i,:)));
tsc(i,kpc) = tsnew(i,kpc);
n_not_nan(i,kpc) = ones(1,length(kpc));
end
n = n_not_nan * n_not_nan';
c = (tsc * tsc') * 1./n;
[lam, pcs, per] = eof(c);
pcs = ones(ntim,1) * (1./(std(pcs))) .* (pcs - ones(ntim, 1)*mean(pcs));
%lds = pcs(:,1:10)'*ts/ntim;
tm = ntim
num = 1;
tem2 = pcs(1:tm,num);
tem2 = (tem2 - mean(tem2)) / std(tem2);
tem = NaN * ones(1, nlat*nlon);
tem(1, kp) = regress_nan(tem2, ts(1:tm, :));
tem = reshape(tem, nlat, nlon);
figure(1)
sd(1)
gcont(tem, [-5:.25:5]);
dc
sd(2)
plot(tem2);
%set(gca, 'XTick', [1:120:1153], 'XTickLabel', [1900:10:2000])
axis([0 (tm+1) -3 3])
%axis([1 ntim+1-500 -3 4])
grid
% Now regress ct* against GR:
g = pcs(:,1);
g = (g - mean(g)) / std(g);
ind = find(~isnan(ctstar(1:500)));
gr = g - corr(g(ind), ctstar(ind)) * ctstar(kptime);
gr = (gr - mean2(gr)) / std(gr(ind));
tem = NaN * ones(1, nlat*nlon);
tem(1,kp) = regress_nan(gr, ts);
tem = reshape(tem, nlat, nlon);
tem = myrunning_ave(myrunning_ave(tem', 2)', 2);
cint = 0.1;
sd(1)
gcont(tem, [-5:cint:5])
dc
title(['COADS: GR pattern'])
xlabel(['Contour Interval: ' num2str(cint) ' K (std)^-^1'])
sd(2)
plot(gr)
set(gca, 'XTick', [1:60:626], 'XTickLabel', [1950:5:2000])
axis([0 ntim+1 -3 4])
xlabel('Year: ticks correspond to January')
ylabel('STD')
title(['Low Pass GR time series'])
grid
% Check out other possibilities:
grlp = myrunning_ave(gr, 25);
grlp = myrunning_ave(grlp, 37);
grlp = (grlp - mean2(grlp)) / std(grlp(find(~isnan(grlp))));
cint = 0.1;
tem = NaN * ones(1, nlat*nlon);
tem(1, kp) = regress_nan(grlp, ts);
tem = reshape(tem, nlat, nlon);
%tem = myrunning_ave(tem, 2);
%tem = myrunning_ave(tem', 2);
%tem = tem';
sd(1)
gcont(tem, [-5:cint:5])
dc
title(['COADS2x2: Low Pass GR pattern'])
xlabel(['Contour Interval: ' num2str(cint) ' K (std)^-^1'])
sd(2)
plot(grlp)
set(gca, 'XTick', [1:60:626], 'XTickLabel', [1950:5:2000])
axis([0 ntim+1 -2 2])
xlabel('Year: ticks correspond to January')
ylabel('STD')
title(['Low Pass GR time series'])
grid
cd /home/disk/hayes2/dvimont/data/Plots
grhp = gr - myrunning_ave(myrunning_ave(gr, 25), 37);
grhp = (grhp - mean(grhp)) / std(grhp);
cint = 0.1;
tem = NaN * ones(1, nlat*nlon);
tem(1, kp) = regress_nan(grhp, ts);
tem = reshape(tem, nlat, nlon);
tem = myrunning_ave(myrunning_ave(tem', 2), 2)';
sd(1)
gcont(tem, [-5:cint:5])
dc
title(['COADS2x2: High Pass GR pattern'])
xlabel(['Contour Interval: ' num2str(cint) ' K (std)^-^1'])
sd(2)
plot(grhp)
set(gca, 'XTick',[49:60:409],'XTickLabel',[65:5:100],'YTick',[-3:3])
axis([0 409 -3 3.5])
xlabel('Year: ticks correspond to January')
ylabel('STD')
title(['High Pass GR time series'])
grid
cd /home/disk/hayes2/dvimont/data/Plots
% Check lag/lead relationship:
corr_grct = [];
for i = -24:12;
corr_i = corr(ctstar, gr, i);
corr_grct = [corr_grct corr_i];
end
subplot(2,1,1)
bar(-24:12, corr_grct)
grid
title(['Lagged Correlations between GR and CT*'])
xlabel(['Months: Negative => CT* leads GR; Positive => CT* lags GR'])
ylabel(['Correlation Coefficient'])
axis([-25 13 -0.2 0.5])
corr_grct = [];
corr_hpgrct = [];
for i = -24:12;
corr_i = corr(ctstar, grhp, i);
corr_grct = [corr_grct corr_i];
end
subplot(2,1,2)
bar(-24:12, corr_grct)
grid
title(['Lagged Correlations between High Pass GR and CT*'])
xlabel(['Months: Negative => CT* leads GR; Positive => CT* lags GR'])
ylabel(['Correlation Coefficient'])
axis([-25 13 -0.5 0.8])
cd /home/disk/hayes2/dvimont/data/Plots
% Just for the heck of it, try rotating the loadings.
nkp = 80; % Allows 50% var explained
lds = zeros(nkp, size(ts, 2));
for i = 1:nkp;
lds(i,:) = regress_nan(pcs(:,i)', ts);
end
wgt = diag(lds * 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';