Global Index (short | long) | Local contents | Local Index (short | long)
ctlim = [110 300 -30 30]; ctlim = [110 300 -72.5 67.5];
This script calls | |
---|---|
clear for index = 1:3 cd /home/disk/hayes2/dvimont/csiro/data if index == 1; tim = [101:550]; elseif index == 2; tim = [551:1000]; else; tim = [101:1000]; end filin = ['temp_A_L1-10.nc']; nc = netcdf(filin, 'nowrite'); depth = nc{'depth'}(:); lat = nc{'latitude'}(:); lon = nc{'longitude'}(:); ctlim = [180 270 -6 6]; [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))); lat = lat(yk); lon = lon(xk); get_global; default_global; FRAME = ctlim; temp = squeeze(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 i = 1:(nlat/2); for j = 1:(nlon/2); yk = 2*(i-1)+[1:2]; xk = 2*(j-1)+[1:2]; temp2(i,j,:) = mean2(mean2(temp(yk,xk,:))); end end temp = shiftdim(temp2, 2); clear temp2; [ntim, nlat, nlon] = size(temp); end temp = reshape(temp, ntim, nlat*nlon); temp = detrend(temp); %temp = temp ./ (ones(ntim, 1) * std(temp)); temp = cosweight(temp, lat); clim = mean(temp); kp = find(~isnan(clim)); nkp = length(kp); temp = temp(:, kp); [lam, per, lds, pcs] = eof_dan(temp); figure(2); figure_orient; get_global; default_global; FRAME = ctlim; cint = 0.1; clev = -1:cint:1; num = 1; tem = NaN*ones(nlat*nlon, 1); tem(kp) = lds(:,num); tem = reshape(tem, nlat, nlon); sd(1) mcont(tem, clev, 'giso') title(['EOF' num2str(num) ' of level ' num2str(num) ... ' TEMP, ' num2str(round(per(num))) '% Variance explained, ' ... 'Years ', num2str(tim(1)), ' to ', num2str(tim(ntim))]); xlabel(['Contour Interval: ' num2str(cint) ' K std^-^1']); sd(2) plot(tim, pcs(:,num)); title(['PC' num2str(num)]) xlabel('Year'); ylabel('STD'); axis([min(tim)-1 max(tim)+1 -3 3]); grid on; cd /home/disk/tao/dvimont/matlab/CSIRO/Plots2 figure(2);% 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') 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(['LEVEL 1 TEMP: PC', num2str(num)]); subplot(4,1,4) xlabel('Years') cd /home/disk/tao/dvimont/matlab/CSIRO/Plots2 index if index == 1; pc101 = pcs; elseif index == 2; pc551 = pcs; end end % Calculate power spectrum for pcs ctstar = detrend(squeeze(mean2(mean2(shiftdim(temp, 1))))); ctstar = ctstar./std(ctstar); num = 1 for i = 1:3; if i == 1; % ctann = pc101(:,num); ctann = ctstar(1:450); tit = ['Years 101 to 550']; elseif i == 2; % ctann = pc551(:,num); ctann = ctstar(451:900); tit = ['Years 551 to 1000']; elseif i == 3; % ctann = pcs(:,num); ctann = ctstar(1:900); tit = ['Years 101 to 1000']; end nx = length(ctann); nfft = 128; 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)); if i ~= 3; rner1 = rn * 2.64; rner5 = rn * 2.01; else rner1 = rn * 2.05; rner5 = rn * 1.68; end dofx = nx / n2; figure(1) subplot(3,1,i) semilogy(f2, p(:,1), 'b-', f2, rn, 'b-', ... f2, rner5, 'b--', f2, rner1, 'b-.') % tim = 2:(n2+1); f2 = f2(tim); p = p(tim,:); rn = rn(tim); % rner1 = rner1(tim); rner5 = rner5(tim); % plot(1*log(f2), f2'.*p(:,1), 'b-', 1*log(f2), f2.*rn, 'b-', ... % 1*log(f2), f2.*rner5, 'b--', 1*log(f2), f2.*rner1, 'b-.') set(gca, 'YTick', [.01 .1 1 10 100]); axis([0 0.5 .05 15]) grid title(['Power Spectrum for PC' num2str(num) ', 30s to 30n, ' tit]); % title(['Power Spectrum for CT, ' tit]); xlabel(['Frequency: yr^-^1; Degrees of Freedom: ' ... num2str(round(dofx)) ';' ... ' Bandwidth: 7.8 x 10^-^3 yr^-^1']) l=legend('PC1 Spectrum', 'Red Noise', '95% Confidence', '99% Confidence'); end cd /home/disk/tao/dvimont/matlab/CSIRO/Plots2