Global Index (short | long) | Local contents | Local Index (short | long)
tim = [551:1000];
| This script calls | |
|---|---|
clear
cd /home/disk/hayes2/dvimont/csiro/data
filin = 'temp_A_L1-10.nc';
tim = [101:550];
nc = netcdf(filin, 'nowrite');
depth = nc{'depth'}(:);
lat = nc{'latitude'}(:);
lon = nc{'longitude'}(:);
ctlim = [110 300 -60 60];
[xk, yk] = keep_var(ctlim, lon, lat);
temp = nc{'temp'}(tim,1,yk,xk);
mv = nc{'temp'}.missing_value(:);
nc = close(nc);
[ntim, nlev, nlat, nlon] = size(temp);
temp(find(temp == mv)) = NaN * ones(size(find(temp == mv)));
lat1 = lat(yk); lon1 = lon(xk);
temp = reshape(temp, ntim, nlat*nlon);
get_global; default_global; FRAME = ctlim;
% Look at location of 20degree isotherm
if 0;
[xk, yk] = keep_var([110 300 -6 6], lon, lat);
tim = [101:550];
tim = [551:1000];
tem = squeeze(mean(shiftdim(mean(temp(500:510,:,yk,xk)), 2)))';
contour(lon, -.01*depth, tem, [0:1:20], 'k')
end
% We'll start with using the 20degree isotherm
%
% First, find min and max latitude of the 20deg isotherm
clim = squeeze(mean(temp(:,:,:,:)));
tiso = 20;
latind = [];
for i = 1:nlon;
tem = [max(find(clim(1,:,i) >= tiso)) min(find(clim(1,:,i) >= tiso))];
latind = [max([latind tem]) min([latind tem])];
end
if rem(length(latind(2):latind(1)), 2);
latind(1) = latind(1)+1;
end
latind = latind(2):latind(1);
temp = temp(:,:,latind,:);
clim = clim(:,latind,:);
[ntim, nlev, nlat, nlon] = size(temp);
if 0;
temp = shiftdim(temp, 2);
clim = shiftdim(clim, 1);
temp2 = NaN * ones(ntim, nlev, nlat/2, nlon/2);
clim2 = NaN * ones(nlev, nlat/2, nlon/2);
for i = 1:nlat/2;
indy = 2*(i-1)+[1:2];
for j = 1:nlon/2;
indx = 2*(j-1)+[1:2];
temp2(:,:,i,j) = squeeze(mean2(mean2(temp(indy,indx,:,:))));
clim2(:,i,j) = squeeze(mean2(mean2(clim(indy,indx,:))));
end
end
lat = lat(latind);
lat = mean([lat(1:2:nlat)' ; lat(2:2:nlat)']);
lon = mean([lon(1:2:nlon)' ; lon(2:2:nlon)']);
else
temp2 = temp;
clim2 = clim;
end
[ntim, nlev, nlat, nlon] = size(temp2);
kp20 = find(clim2(1,:) >= tiso);
clim2 = reshape(clim2, nlev, nlat*nlon);
temp2 = reshape(temp2, ntim, nlev, nlat*nlon);
clim2 = clim2(:,kp20);
temp2 = temp2(:,:,kp20);
[ntim, nlev, nkp] = size(temp2);
ind = [];
for j = 1:nkp;
if isempty(find(diff(clim2(:,j)) > 0));
ind = [ind j];
end
end
temp20 = NaN * ones(ntim, nkp);
for i = 1:ntim;
for j = ind;
temp20(i,j) = interp1(clim2(:,j), squeeze(temp2(i,:,j)), 20);
end
end
% cd /home/disk/hayes2/dvimont/csiro/matlab_data
% save 20degiso_yr551-1000_temp_noave.mat temp20 kp20 lat lon tim
biff = 3;
ind = 1;
cd /home/disk/hayes2/dvimont/csiro/matlab_data
if ind == 1;
load 20degiso_yr101-550_temp.mat
lag = -60;
ttit = ['101 - 550'];
elseif ind == 2;
load 20degiso_yr551-1000_temp_noave.mat
lag = 90;
ttit = ['551 - 1000'];
end
ctlim = [110 300 -75 65];
ntim = length(tim); nlat = length(lat); nlon = length(lon);
temp3 = NaN * ones(ntim, nlat*nlon);
temp3(:,kp20) = temp20;
temp3 = detrend(temp3);
temp3 = cosweight(temp3, lat);
kp = find(~isnan(mean(temp3)));
temp3 = temp3(:,kp);
% [lam1, per1, lds1, pcs1] = eof_dan(temp3);
% [lam2, per2, lds2, pcs2] = eof_dan(temp3);
if biff == 1;
yr1 = 10; yr2 = 10; tit = ['HP ( < 10 yr ) '];
elseif biff == 2;
yr1 = 20; yr2 = 20; tit = ['HP ( < 20 yr ) '];
elseif biff == 3;
yr1 = 10; yr2 = 60; tit = ['LP ( > 10 yr ) '];
elseif biff == 4;
yr1 = 20; yr2 = 60; tit = ['LP ( > 20 yr ) '];
elseif biff == 5;
yr1 = 5; yr2 = 20; tit = ['HP ( 8 yr ) '];
elseif biff == 6;
yr1 = 10; yr2 = 50; tit = ['BP ( 10-50 yr ) '];
elseif biff == 7;
tit = ['Unfiltered Data '];
end
[b1, a1] = butter(6, 2/yr1); [b2, a2] = butter(6, 2/yr2);
if biff < 3 | biff == 5;
temp3 = temp3 - filtfilt(b1, a1, temp3); disp('HP Filtered');
elseif biff > 2 & biff < 5;
temp3 = filtfilt(b1, a1, temp3); disp('LP Filtered');
elseif biff ~= 7
temp3 = filtfilt(b1, a1, temp3) - filtfilt(b2, a2, temp3); disp('BP Filtered');
end
% Do complex eof:
[lam, lds, pcs, per] = complex_eof(temp3);
j = sqrt(-1);
npt = length(kp);
pcs = sqrt(2)*pcs ./ (ones(ntim,1) * std(pcs));
% Look at different phases 11111111111111111111111111111111111111111111111111111
temp3 = NaN * ones(ntim, nlat*nlon);
temp3(:,kp20) = temp20;
%temp2 = temp; XAX = lon1; YAX = lat1;
temp2 = temp3; XAX = lon; YAX = lat;
%temp2 = -1*hflx; XAX = lon2; YAX = lat2;
nlat = length(YAX); nlon = length(XAX);
temp2 = detrend(temp2);
if biff < 3 | biff == 5;
temp2 = temp2 - filtfilt(b1, a1, temp2);% - filtfilt(b2, a2, temp2);
elseif biff > 3 & biff < 5;
temp2 = filtfilt(b1, a1, temp2);
elseif biff == 6;
temp2 = filtfilt(b1, a1, temp2) - filtfilt(b2, a2, temp2);
end
stdtemp = std(temp2);
dofx = floor(ntim/10 - 3);
FRAME = [ctlim(1:2) floor(min(YAX)) ceil(max(YAX))];
figure(3); clf; figure_orient;
lag = 0; lg = lag*pi/180; lg2 = 1;
num = 1; lind = 1;
nfrm = 6;
cint = 0.05; clev = [-9:cint:9];
for i = 1:nfrm;
wgt = conj(exp(j * ((i-1) * pi/(lg2*nfrm) + lg) ));
temtim = squeeze(real(wgt .* pcs(:,num)));
tem = temtim' * temp2 ./ ntim;
% ccoef = tem./stdtemp;
% tstat = ccoef ./ (sqrt((1-ccoef.^2)/dofx));
% score = tscore(dofx, 0.5);
tem = reshape(tem, nlat, nlon);
% tstat = reshape(tstat, nlat, nlon);
subplot(3,2,i);
cla;
% gshade2(abs(tstat), score); hold on;
gcont(tem, clev); hold off;
axis([ctlim(1:2) -60 60]);
dc
title(['Phase = ' num2str(((i-1)*180)/(lg2*nfrm)+lag)]);
xlabel(['Contour Interval: ' num2str(cint) ' K std^-^1'])
end
subplot(3,2,3);
ylabel([tit ': < Surface Temp, CPC1 of 20 deg. T''>; Years ' ttit]);
ylabel([tit ': < 20 deg. T'', CPC1 of 20 deg. T'' >; Years ' ttit]);
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots5/yr101-550
tit1 = tit;
% print -dps2 T_L7_20deg_CEOF1_yr551-1000.ps
% Plot the CPC's 2222222222222222222222222222222222222222222222222222222222222
if ind == 1; tim = 101:550; else; tim = 551:1000; end;
figure(2); figure_orient;
for m = 1:3
subplot(4,1,m)
tind = ((ntim/3)*(m-1)+[1:(ntim/3)]);
plot(tim(tind), real(pcs(tind, num)), '-k', ...
tim(tind), imag(pcs(tind, num)), '--k');
axis([min(tim(tind)-1) max(tim(tind)+1) -3 3]);
set(gca, 'YTick', [-3:3.5], 'XTick', [(tim(1)-1):10:max(tim)]);
grid
end
subplot(4,1,1)
title([tit ': CPC1 20 DEG ISOTHERM: ' num2str(round(per(1)))...
'% Variance Explained']);
subplot(4,1,2)
ylabel(['Solid = REAL, Dashed = IMAG']);
subplot(4,1,3)
xlabel('Years')
if 1;
subplot(4,2,7);
plot(real(pcs(:,1)), imag(pcs(:,1)), '.k');
axis square
axis([-4 4 -4 4]); grid on
xlabel('Re(CPC1)'); ylabel('Im(CPC1)');
end;
subplot(4,2,7);
binnum = 24;
binind = [(-1*(pi-(pi/binnum))):(2*pi/binnum):(pi-(pi/binnum))];
ph = atan2(imag(pcs(:,1)), real(pcs(:,1)));
hgram = hist(ph, binind);
bar((180/pi)*binind, hgram);
axis([-180 180 0 30])
set(gca, 'XTick', -180:90:180);
grid on
xlabel('PHASE (deg)');
ylabel('COUNT');
subplot(4,2,8); nind = 1;
tim = 1:450;
v1 = real(pcs(tim,1)); v2 = imag(pcs(tim,1));
lab = ['Re(CPC1)'; 'Im(CPC1)'];
tind = (-8*nind):nind:(8*nind);
a = zeros(length(tind),1);
for lag = tind;
a((1/nind)*(lag+max(tind))+1) = corr(v1, v2, lag);
end
bar(tind, a);
title(['< ' lab(1,:) ', ' lab(2,:) ' >'])
% xlabel([lab(1,:) ' leads ' lab(2,:) ' | ' lab(2,:) ' leads ' lab(1,:)])
xlabel([lab(1,:) ' leads | ' lab(2,:) ' leads '])
axis([min(tind-1) max(tind+1) -1.1 1.1])
set(gca, 'XTick', min(tind):2*nind:max(tind), 'YTick', -1:.25:1)
grid on
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots5
% Look at real and imaginary components
default_global;
FRAME = [ctlim(1:2) floor(min(lat)) ceil(max(lat))];
figure(2); figure_orient;
num = 1;
subplot(3,1,1)
tem = NaN * ones(1, nlat*nlon);
tem(kp) = real(lds(:,num));
tem = reshape(tem, nlat, nlon);
cint = 0.1; clev = [-1:cint:1];
gcont(tem, clev); dc;
subplot(3,1,2)
tem = NaN * ones(1, nlat*nlon);
tem(kp) = imag(lds(:,num));
tem = reshape(tem, nlat, nlon);
cint = 0.1; clev = [-1:cint:1];
gcont(tem, clev); dc;
subplot(6,1,5);
tind = 1:ntim/2;
plot(tim(tind), real(pcs(tind,num)), '-k', ...
tim(tind), imag(pcs(tind,num)), '--k')
axis([min(tim(tind)-1) max(tim(tind)+1) -3 3]);
grid on;
set(gca, 'XTick', [(tim(min(tind))-1):25:(tim(max(tind)))]);
subplot(6,1,6);
tind = [1:ntim/2] + ntim/2;
plot(tim(tind), real(pcs(tind,num)), '-k', ...
tim(tind), imag(pcs(tind,num)), '--k')
axis([min(tim(tind)-1) max(tim(tind)+1) -3 3]);
grid on;
set(gca, 'XTick', [(tim(min(tind))-1):25:(tim(max(tind)))]);
% For fun, a movie.
figure(1); clf;
ph = atan2(imag(pcs(:,1)), real(pcs(:,1)));
ph2 = interp(ph(18:29), 3);
ph2 = 0:pi/12:2*pi;
nmovie = length(ph2);
M = moviein(nmovie);
for i = 1:nmovie;
wgt = conj(exp(j * ((i-1) * 2 * pi/nmovie + lg) ));
temtim = squeeze(real(wgt .* pcs(:,num)));
temtim = (temtim - mean(temtim)) ./ std(temtim);
tem = temtim' * temp3 ./ ntim;
tem = reshape(tem, nlat, nlon);
surf(XAX,YAX,tem);% dc;
axis([110 300 -40 40 -.4 .4]);
view([217.5 40])
M(:,i) = getframe;
end
movie(M, 4, 8)
% Look at EOF
[lam, per, lds, pcs] = eof_dan(temp3);
default_global;
FRAME = [ctlim(1:2) floor(min(lat)) ceil(max(lat))];
figure(1); figure_orient;
num = 1;
sd(1);
tem = NaN * ones(1, nlat*nlon);
tem(kp) = real(lds(:,num));
tem = reshape(tem, nlat, nlon);
cint = 0.1; clev = [-1:cint:1];
gcont(tem, clev); dc;
title(['EOF1 of T'' along Clim. 20 deg. Isotherm; ' ...
num2str(round(per(1))) '% Variance Explained; yr ' ttit]);
xlabel(['Contour Interval: ' num2str(cint) ' K std^-^1']);
sd(2);
plot(tim, pcs(:,1), 'k');
grid on;
axis([min(tim-1) max(tim+1) -3 3]);
title(['PC1']);
ylabel('STD');
xlabel('Years');
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots5
% Look at heat flux
default_global;
FRAME = [110 300 floor(min(lat)) ceil(max(lat))];
cd /home/disk/hayes2/dvimont/csiro/data
nc = netcdf('heat_flux_A_1000_years.nc', 'nowrite');
lat = nc{'latitude'}(:);
lon = nc{'longitude'}(:);
[xk, yk] = keep_var(FRAME, lon, lat);
temp3 = nc{'stfht'}(tim, yk, xk);
mv = nc{'stfht'}.missing_value(:);
nc = close(nc);
lon = lon(xk); lat = lat(yk);
temp3(find(temp3 == mv)) = NaN * ones(size(find(temp3 == mv)));
[ntim, nlat, nlon] = size(temp3);
temp3 = reshape(temp3, ntim, nlat*nlon);
temp3 = detrend(temp3);
% Try the eof on the entire data set
cd /home/disk/hayes2/dvimont/csiro/matlab_data
load 20degiso_yr101-550_temp_noave.mat
load 20degiso_yr551-1000_temp_noave.mat
ntim = length(tim); nlat = length(lat); nlon = length(lon);
temp3 = NaN * ones(ntim, nlat*nlon);
temp3(:,kp20) = temp20;
temp3 = detrend(temp3);
temp3 = cosweight(temp3, lat);
kp = find(~isnan(mean(temp3)));
temp3 = temp3(:,kp);
temp1 = temp3; kp1 = kp;
[lam1, per1, lds1, pcs1] = eof_dan(temp1);
temp2 = temp3; kp2 = kp;
[lam2, per2, lds2, pcs2] = eof_dan(temp2);
kp = kp1(find(ismember(kp1, kp2)));
temp3 = [temp1(:, find(ismember(kp1, kp2))); temp2(:, find(ismember(kp2, kp1)))];
[lam3, per3, lds3, pcs3] = eof_dan(temp3);