Global Index (short | long) | Local contents | Local Index (short | long)
ctlim = [-0.1 17.5 47.5 80]; ctlim = [300 360 47.5 80]; ctlim = [110 300 -30 30];
This script calls | |
---|---|
clear cd /home/disk/hayes2/dvimont/csiro/data tim = [101:550]; filin = ['temp_A_L1-10.nc']; nc = netcdf(filin, 'nowrite'); depth = nc{'depth'}(:); lat = nc{'latitude'}(:); lon = nc{'longitude'}(:); ctlim = [110 300 -67.5 60]; % ctlim = [110 300 12.5 62.5]; % ctlim = [110 300 -22.5 22.5]; % ctlim = [110 300 22.5 62.5]; [xk, yk] = keep_var(ctlim, lon, lat); temp = nc{'temp'}(tim,:,yk,xk); mv = nc{'temp'}.missing_value(:); nc = close(nc); lev2 = 1:3:10; temp = temp(:,lev2,:,:); temp(find(temp == mv)) = NaN * ones(size(find(temp == mv))); lat = lat(yk); lon = lon(xk); [ntim, nlev, nlat, nlon] = size(temp); % Average data if necessary: if 1 temp = shiftdim(temp, 1); lat = mean([lat(1:2:nlat)'; lat(2:2:nlat)'])'; lon = mean([lon(1:2:nlon)'; lon(2:2:nlon)'])'; for l = 1:nlev; for i = 1:(nlat/2); % for i = 1:nlat; for j = 1:(nlon/2); yk = 2*(i-1)+[1:2]; xk = 2*(j-1)+[1:2]; % yk = i; xk = 2*(j-1)+[1:2]; temp2(l,i,j,:) = mean2(mean2(squeeze(temp(l,yk,xk,:)))); % temp2(l,i,j,:) = mean2(squeeze(temp(l,yk,xk,:))); end end end end temp = shiftdim(temp2, 3); clear temp2; [ntim, nlev, nlat, nlon] = size(temp); % Do this on the correlation matrix: temp = reshape(temp, ntim, nlev*nlat*nlon); for i = 1:nlev*nlat*nlon; temp(:,i) = temp(:,i)./std(temp(:,i)); end temp = reshape(temp, ntim, nlev, nlat, nlon); % Weight by depth if 0; middepth = [0; (depth([lev2])+depth(lev2+1)) / 2]; wgt = diff(middepth)/100; wgt = wgt./mean(wgt); temp = shiftdim(temp, 1); temp = reshape(temp, nlev, nlat*nlon*ntim); temp = temp .* (wgt * ones(1, nlat*nlon*ntim)); temp = reshape(temp, nlev, nlat, nlon, ntim); temp = shiftdim(temp, 3); end for i = 1:nlev; temp(:,i,:,:) = cosweight(squeeze(temp(:,i,:,:)), lat); end; temp = reshape(temp, ntim, nlev*nlat*nlon); temp = detrend(temp); clim = mean(temp); kp = find(~isnan(clim)); npt = length(kp); temp = temp(:, kp); % Lowpass filter if 1 yr = 8; [b,a]=butter(6,(2/yr)); temp = temp - filtfilt(b, a, temp); end if 0; c = temp * temp' / npt; [lam, pcs, per] = eof(c); pc10 = pcs(:,1:10); pc10 = (pc10 - ones(ntim, 1)*mean(pc10)) ./ (ones(ntim, 1)*std(pc10)); ld10 = pc10' * temp ./ ntim; cd /home/disk/hayes2/dvimont/csiro/matlab_data %save temp_3deof_pac.mat pc10 ld10 kp lat lon lev2 lam per clear c pcs end cd /home/disk/hayes2/dvimont/csiro/matlab_data [lam, lds, pcs, per] = complex_eof(temp); %save 6yrLP_PAC_NATL_CEOF.ps lam lds pcs per kp lat lon lev2 depth kp kp2 kp3 lat lon lat2 lon2 lat3 lon3 %save LPtemp_3dceof_22.5to62.5pac.mat lam lds pcs per kp lat lon lev2 depth %save HPtemp_3dceof_12.5to62.5pac.mat lam lds pcs per kp lat lon lev2 depth %save 8yrHPtemp_3dceof_-67.5to60pac.mat lam lds pcs per kp lat lon lev2 depth load tem_latlon.mat % % % clear cd /home/disk/hayes2/dvimont/csiro/matlab_data load temp_3dceof_pac.mat get_global; default_global; FRAME = [110 300 -75 65]; ntim = size(pcs, 1); npt = length(kp); nlev = length(lev2); nlat = length(lat); nlon = length(lon); % Weight the CPCs and CEOFs accordingly weight1 = std(real(pcs)); weight2 = std(imag(pcs)); pcs = real(pcs).*(ones(ntim, 1)*(1./weight1)) + ... j*imag(pcs).*(ones(ntim, 1)*(1./weight2)); lds = real(lds).*(ones(npt, 1)*weight1) + ... j*imag(lds).*(ones(npt, 1)*weight2); % Reshape temperature field temp = NaN*ones(10, nlev*nlat*nlon); temp(:,kp) = lds'; temp = reshape(temp, 10, nlev, nlat, nlon); % Unweight the depth dimension middepth = [0; (depth([lev2])+depth(lev2+1)) / 2]; wgt = diff(middepth)/100; wgt = wgt./mean(wgt); temp = shiftdim(temp, 1); temp = reshape(temp, nlev, nlat*nlon*10); temp = temp ./ (wgt * ones(1, nlat*nlon*10)); temp = reshape(temp, nlev, nlat, nlon, 10); temp = shiftdim(temp, 3); for lev = 1:5 figure(lev); figure_orient; num = 1; tem = real(squeeze(temp(num, lev, :, :))); cint = 0.05; clev = [-.5:cint:.5]; sp(1) if ismap(gca); clma; end; mcont(tem, clev, 'giso') title(['3D TEMP: Real CEOF' num2str(num) ', ' ... num2str(depth(lev2(lev))/100) 'm DEPTH, ' ... num2str(round(per(num))) '% Variance Explained']); xlabel(['Contour Interval: ' num2str(cint) ' K std^-^1']); tem = imag(squeeze(temp(num, lev, :, :))); cint = 0.05; clev = [-.5:cint:.5]; sp(2) if ismap(gca); clma; end; mcont(tem, clev, 'giso') title(['3D TEMP: Imag CEOF' num2str(num) ', ' ... num2str(depth(lev2(lev))/100) 'm DEPTH, ' ... num2str(round(per(num))) '% Variance Explained']); xlabel(['Contour Interval: ' num2str(cint) ' K std^-^1']); end for lev = 5:-1:1; figure(lev); end; figure(6);% figure_orient; for m = 1:4 subplot(4,1,m) tind = ((ntim/4)*(m-1)+[1:(ntim/4)]); plot(tind, real(pcs(tind, num)), '-k', ... tind, imag(pcs(tind, num)), '--k'); axis([min(tind-1) max(tind+1) -3 3]); set(gca, 'YTick', [-3:3], 'XTick', [0:25:ntim]); grid end subplot(4,1,1) title(['3D TEMP: CPC', num2str(num), ... ': Dashed = IMAG, Solid = REAL']); subplot(4,1,4) xlabel('Years') cd /home/disk/tao/dvimont/matlab/CSIRO/Plots % Out of curiosity, look at spectrum of CPC: for i = 1:2; if i == 1; ctann = real(pcs(1:450,num)); elseif i == 2; ctann = imag(pcs(1:450,num)); end nx = length(ctann); nfft = 100; noverlap = nfft/2; [p, f] = spectrum(ctann, nfft, noverlap); n2 = nfft/2; f2 = 0.5 * (0:n2)/n2; rho = (corr(ctann, ctann, 1) + sqrt(corr(ctann, ctann, 2))) / 2; % + ... % (corr(ctann, ctann, 3)^(1/3))) / 3; rn = (1-rho^2) ./ (1-2*rho*cos((pi/n2)*(0:n2))+rho^2); rn = rn / mean(rn); rn = rn * mean(p(:,1)); rner1 = rn * 2.05; rner5 = rn * 1.68; dofx = nx / n2; figure(1) sp(i) semilogy(f2, p(:,1), 'b-', f2, rn, 'b-', ... f2, rner5, 'b--', f2, rner1, 'b-.') set(gca, 'YTick', [.01 .1 1 10 100]); axis([0 0.5 .05 15]) grid title(['Power Spectrum for CPC' num2str(num) ', 3D TEMP']); xlabel(['Frequency: yr^-^1; Degrees of Freedom: ' ... num2str(round(dofx)) ';'... ' Bandwidth: 7.8 x 10^-^3 yr^-^1']) l=legend('CPC Spectrum', 'Red Noise', '95% Confidence', '99% Confidence'); end