Global Index (short | long) | Local contents | Local Index (short | long)
load ML_SLP_eof_shem_-60to-20.mat load ML_SLP_eof_spac.mat;
| This script calls | |
|---|---|
clear
dofilt = 0;
data
load ML_SLP_eof_npac.mat;
back
mrpcs = rpcs; mlpcs = lpcs; mhpcs = hpcs;
mrper = rper; mlper = lper; mhper = hper;
mrlds = rlds; mllds = llds; mhlds = hlds;
mrlam = rlam; mllam = llam; mhlam = hlam;
if dofilt
uselds = 1*mllds; usepcs = 1*mlpcs;
% uselds = 1*mhlds; usepcs = 1*mhpcs;
uselims = lims;
else
uselds = 1*mrlds; usepcs = 1*mrpcs;
uselims = lims;
end
% Get ppcs
tim = 101:1000;
slpc = getflx('psl', lims, tim);
csirod
load ML_ANN_slp.mat;
[xk, yk] = keep_var(uselims, lon, lat);
slp = slp(:, yk, xk);
back
if dofilt;
[b, a] = butter(9, 2/9);
slpc = filtfilt(b, a, detrend(slpc));
slp = filtfilt(b, a, detrend(slp));
% [b, a] = butter(9, 2/10);
% slpc = detrend(slpc)-filtfilt(b, a, detrend(slpc));
% slp = detrend(slp)-filtfilt(b, a, detrend(slp));
else
slpc = detrend(slpc);
slp = detrend(slp);
end;
% Normalize uselds;
wgt = sqrt(diag(uselds'*uselds));
uselds = uselds ./ (ones(size(uselds, 1), 1)*wgt');
% Project lds onto slp
[ntim, nlat, nlon] = size(slpc);
ppcs = reshape(slpc, ntim, nlat*nlon) * uselds;
[ntim, nlat, nlon] = size(slp);
ml_ppcs = reshape(slp, ntim, nlat*nlon) * uselds;
lims = [-0.1 360 -90 90];
sstc = getnc('temp', lims, 1, tim);
slpc = getflx('psl', lims, tim);
%txc = getnc('z250', lims, 1, tim);
txc = getnc('taux', lims, 1, tim);
csirod
load ML_ANN_slp.mat;
load ML_ANN_sst.mat; sst = sst2; clear sst2;
%load z250_ML_annave.mat; txann = z250;
load taux_ML_annave.mat
back
lims = [-0.1 360 -90 90];
[xk, yk] = keep_var(lims, lon, lat);
slp = slp(:, yk, xk);
sst = sst(:, yk, xk);
txann = txann(:, yk, xk);
ctc = getct;
[xkc, ykc] = keep_var([180 270 -6 6], lon, lat);
ctm = squeeze(mean2(mean2(shiftdim(sst(:,ykc,xkc), 1))));
ctc = detrend(ctc); ctm = detrend(ctm);
if dofilt;
ctm = filtfilt(b, a, ctm);
ctc = filtfilt(b, a, ctc);
% ctm = ctm - filtfilt(b, a, ctm);
% ctc = ctc - filtfilt(b, a, ctc);
end
% Get regression maps
if dofilt;
slp = filtfilt(b, a, detrend(slp)); slpc = filtfilt(b, a, detrend(slpc));
sst = filtfilt(b, a, detrend(sst)); sstc = filtfilt(b, a, detrend(sstc));
txann = filtfilt(b, a, detrend(txann)); txc = filtfilt(b, a, detrend(txc));
% slp = detrend(slp-filtfilt(b, a, slp)); slpc = detrend(slpc-filtfilt(b, a, slpc));
% sst = detrend(sst-filtfilt(b, a, sst)); sstc = detrend(sstc-filtfilt(b, a, sstc));
% txann = detrend(txann-filtfilt(b, a, txann)); txc = detrend(txc-filtfilt(b, a, txc));
else
slp = detrend(slp); slpc = detrend(slpc);
sst = detrend(sst); sstc = detrend(sstc);
txann = detrend(txann); txc = detrend(txc);
end
num = 2; pm = sign(corr(usepcs(:,num), ctm));
[reg1m, c1m] = regress_eof(slp, pm*usepcs(:,num), 0);
[reg2m, c2m] = regress_eof(sst, pm*usepcs(:,num), 0);
[reg3m, c3m] = regress_eof(txann, pm*usepcs(:,num), 0);
[reg1c, c1c] = regress_eof(slpc, pm*ppcs(:,num), 0);
[reg2c, c2c] = regress_eof(sstc, pm*ppcs(:,num), 0);
[reg3c, c3c] = regress_eof(txc, pm*ppcs(:,num), 0);
% Plot the data
[lat, lon, depth, lm] = getll('temp', lims);
default_global; FRAME = [105 295 -60 60];
figure(2*dofilt+1); fo(1); clf;
sptalk(4,2,1);
gcont(reg1m, 0.2);
dc2(lm, 0.5, -1000);
color_shade(squeeze(c1m.^2), 0.25, 0.7*[1 1 1]);
grid on;
title(['ML: Regression on LP9\_PC' num2str(num)]);
ylabel('SLP');
sptalk(4,2,2);
gcont(reg1c, 0.2);
dc2(lm, 0.5, -1000);
color_shade(squeeze(c1c.^2), 0.25, 0.7*[1 1 1]);;
grid on;
title(['COUP: Regression on LP9\_PPC' num2str(num)]);
sptalk(4,2,3); FRAME = [105 295 -60 60];
gcont(reg2m, 0.05);
dc2(lm, 0.5, 100);
color_shade(squeeze(c2m.^2), 0.25, 0.7*[1 1 1]);
ylabel('SST');
sptalk(4,2,4);
gcont(reg2c, 0.05);
dc2(lm, 0.5, 100);
color_shade(squeeze(c2c.^2), 0.25, 0.7*[1 1 1]);
[lat1, lon1] = getll('taux', lims);
sptalk(4,2,5); FRAME = [105 295 -60 60];
gcont(10*reg3m, 0.02); % use .02 for wind stress
dc2(lm, 0.5, 100);
color_shade(squeeze(c3m.^2), 0.25, 0.7*[1 1 1]);
ylabel('TAUX');
sptalk(4,2,6);
XAX = lon1; YAX = lat1;
gcont(reg3c, 0.02);
XAX = lon; YAX = lat;
dc2(lm, 0.5, -1000);
color_shade(squeeze(c3c.^2), 0.25, 0.7*[1 1 1]);
cd ~/Thesis/Chap5
% print -dps2 RAW_ML_COUP_PC_PPC_reg.ps
% print -dps2 LP9_ML_COUP_PC_PPC_reg.ps
%data;
%load RAW_detrend_L1-7_EOF_yr101-1000.mat;
%back
tind = 1:450;
tind = 451:900;
tind = 1:900;
clear a ctmc ctcc
lags = -13:13;
for i = 1:length(lags);
a(i,1) = corr(pm*usepcs(:,num), ctm, lags(i));
a(i,2) = corr(pm*ppcs(tind,num), ctc(tind), lags(i));
ctmc(i) = corr(ctm, ctm, lags(i));
ctcc(i) = corr(ctc, ctc, lags(i));
end
figure(2*dofilt+2); fo(1); clf;
sptalk(5,2,1);
h = bar(lags, a(:,1));
axis([-12 12 -0.35 0.85]);
% axis([-6 6 -0.15 0.45]);
grid on;
set(gca, 'XTick', -10:5:10, 'YTick', -1:.25:1);
% set(gca, 'XTick', -10:2:10, 'YTick', -0.5:.1:1);
xlabel('PC1 leads CT | CT leads PC1');
% ylabel('Lagged Correlation');
title('LP9: PC1 and CT Correlation');
% title('RAW: PC1 and CT Correlation');
set(h, 'facecolor', [.7 .7 .7]);
sptalk(5,2,2);
h = bar(lags, a(:,2));
axis([-12 12 -0.35 0.85]);
% axis([-6 6 -0.15 0.45]);
grid on;
set(gca, 'XTick', -10:5:10, 'YTick', -1:.25:1);
% set(gca, 'XTick', -10:2:10, 'YTick', -0.5:.1:1);
% xlabel('PPC1 leads CT | CT leads PPC1');
title('LP9: PPC1 and CT Correlation');
% title('RAW: PPC1 and CT Correlation');
set(h, 'facecolor', [.7 .7 .7]);
% Do fishz stat on this
lind = find(lags == 0);
rho = a(lind,2)*ctcc;
sigz = 1/sqrt(0.5*(get_dof(ctc)+get_dof(ppcs(:,num)))-3);
z = 0.5*log((1+a(:,2))./(1-a(:,2)));
muz = 0.5*log((1+rho)./(1-rho));
sptalk(5,2,4)
h = bar(lags, z);
hold on
plot(lags, muz, '-k', lags, muz+2.0*sigz, '--k', ...
lags, muz-2.0*sigz, '--k');
hold off
axis([-12 12 -0.6 1.1])
% axis([-6 6 -0.15 .55])
set(gca, 'XTick', -10:5:10, 'YTick', -1:.25:1);
% set(gca, 'XTick', -10:2:10, 'YTick', -0.5:.1:1);
grid on
set(h, 'facecolor', [.7 .7 .7]);
title('Fisher''s Z-statistic for PPC1 and CT');
xlabel('PPC1 leads CT | CT leads PPC1');
pos = get(gca, 'Position');
pos(2) = pos(2) - 0.02;
set(gca, 'Position', pos);
cd ~/Thesis/Chap5
% print -dps2 LP9_PC1_PPC1_CT_corr.ps
% print -dps2 RAW_PC1_PPC1_CT_corr.ps
data;
load HP10_detrend_L1-7_EOF_yr101-1000.mat; hcpcs = -pcs;
load COUP_SLP_eof_nhem.mat; slppcs = hpcs;
back
num = 2;
lags = -10:10;
for i = 1:length(lags);
a(i,1) = corr(ppcs(:,num), ctc, lags(i));
a(i,2) = corr(ppcs(:,num), slppcs(:,num), lags(i));
a(i,3) = corr(ppcs(:,num), hcpcs(:,num), lags(i));
end
figure(2); fo(1);% clf
sptalk(3,2,2);
bar(lags, -a(:,1));
axis([-11 11 -0.5 1]);
grid on;
set(gca, 'XTick', -10:5:10, 'YTick', -0.5:.25:1);
xlabel('PPC1 leads CT | CT leads PPC1');
title('Unfiltered');
sptalk(3,2,4);
bar(lags, a(:,2));
axis([-11 11 -0.5 1]);
grid on;
set(gca, 'XTick', -10:5:10, 'YTick', -0.5:.25:1);
xlabel('PPC1 leads SLPPC1 | SLPPC1 leads PPC1');
sptalk(3,2,6);
bar(lags, -a(:,3));
axis([-11 11 -0.5 1]);
grid on;
set(gca, 'XTick', -10:5:10, 'YTick', -0.5:.25:1);
xlabel('PPC1 leads LPPC1 | LPPC1 leads PPC1');
% Look at power spectra
num = 1;
time1 = ppcs(1:900,num)/1;
time2 = ml_ppcs(1:330,num)/1;
nfft = 64; noverlap = 3*nfft/4;
[p1, f1] = spectrum(time1, nfft, noverlap);
[p2, f2] = spectrum(time2, nfft, noverlap);
f = f1./2;
dof2 = 2*length(time2(:,1)) ./ nfft;
dof1 = 2*length(time1(:,1)) ./ nfft;
figure(3); fo(1); clf;
sptalk(3,1,1);
if i == 1;
h = semilogy(f, p1(:,1), 'o-b', f, p2(:,1), 'v-r', ...
f, 4.35*p2(:,2), '--r', f, p2(:,1)/2.98, '--r');
else
h = semilogy(f, p1(:,1), 'o-b', f, p2(:,1), 'v-r', ...
f, 4.35*p2(:,2), '--r', f, p2(:,1)/2.98, '--r');
end
set(h(1:2), 'linewidth', 2)
axis([f(3) 0.5 30 1000])
set(gca, 'XTick', 0:.05:.5);%, 'YTick', [60:10:100 200:100:1000])
grid on;
title('Power Spectra of NP SLP PC1 and PPC1');
xlabel('Frequency (yrs^-^1)')
legend(h, 'PPC1 (COUP)', 'PC1 (ML)');
cd ~/Thesis/Chap5
print -dpsc2 PC1_PPC1_spectrum.ps
ct = getct; ct = detrend(ct);
csirod
load ML_ANN_sst.mat;
back
[xk, yk] = keep_var([180 360 -6 6], lon, lat);
ct2 = sst2(:,yk,xk);
ct2 = detrend(squeeze(mean2(mean2(shiftdim(ct2, 1)))));
num = 1;
time1 = ct;
time2 = ct2;
nfft = 64; noverlap = 3*nfft/4;
[p1, f1] = spectrum(time1, nfft, noverlap);
[p2, f2] = spectrum(time2, nfft, noverlap);
f = f1./2;
dof2 = 2*length(time2(:,1)) ./ nfft;
dof1 = 2*length(time1(:,1)) ./ nfft;
figure(2); fo(1); clf;
sptalk(3,1,1);
h = semilogy(f, p1(:,1), 'o-b', f, p2(:,1), 'v-r', ...
f, 4.35*p2(:,2), '--r', f, p2(:,1)/2.98, '--r');
set(h(1:2), 'linewidth', 2)
axis([f(3) 0.5 .001 .5])
set(gca, 'XTick', 0:.05:.5)
grid on;
title('Power Spectra of CT index');
xlabel('Frequency (yrs^-^1)')
legend(h, 'CT (COUP)', 'CT (ML)');
cd ~/Thesis/Chap5
print -dpsc CT_ML_COUP_spectrum.ps
% Look at coherence
num = 1;
for biff = 1:1;
if biff == 1;
tim = [1:900];
else
tim = 450*(biff-2)+[1:450];
end
time1 = ppcs(tim,num);
time2 = ct(tim); lsty = '-b';
time3 = ml_ppcs(1:330,num);
time4 = ct2; lsty = '-r';
end
nfft = 64; noverlap = 3*nfft/4;
[p1, f1] = spectrum(time1, time2, nfft, noverlap);
[p2, f1] = spectrum(time3, time4, nfft, noverlap);
f = f1./2;
amp1 = abs(p1(:,3));
amp2 = abs(p2(:,3));
phs1 = atan2(imag(p1(:,3)), real(p1(:,3)));
phs2 = atan2(imag(p2(:,3)), real(p2(:,3)));
figure(3); fo(1); clf;
sptalk(3,1,1);
h=plot(f, p1(:,5), 'o-b', f, p2(:,5), 'v-r');
set(h, 'linewidth', 2);
hold on;
hline(0.28, '-.r');
hline(0.105, '--b');
axis([f(3) 0.5 -0.1 1.1]);
set(gca, 'XTick', 0:.05:.5);
grid on
title('Squared Coherence between PC1 (PPC1) and CT index')
legend(h, 'COUP', 'ML');
sptalk(3,1,2);
h=plot(f, 180/pi*phs1, 'ob', f, 180/pi*phs2, 'vr');
set(h, 'linewidth', 2);
axis([f(3) 0.5 -180 180]);
set(gca, 'XTick', 0:.05:.5);
set(gca, 'YTick', -180:90:180);
grid on
title('Phase lag between PC1 (PPC1) and CT index');
xlabel('Frequency (yrs^-^1)');
legend(h, 'COUP', 'ML');
cd ~/Thesis/Chap5
%print -dpsc2 PC1_PPC1_CT_cospectrum.ps
%sptalk(3,1,2);
%h=plot(f, amp1, '-b', f, amp2, '-r');
%set(h, 'linewidth', 2);
%axis([f(3) 0.5 -1 8]);
%set(gca, 'XTick', 0:.05:.5);
%grid on