Global Index (short | long) | Local contents | Local Index (short | long)
filin = 'histape.cdf.36';
This script calls | |
---|---|
cd /home/disk/hayes2/werner/matlab get_global filin = '36yearly.cdf'; nc = netcdf(filin, 'nowrite'); ts = nc{'TS'}(:); lat = nc{'lat'}(:); lon = nc{'lon'}(:); oro = nc{'ORO'}(:); nc = close(nc); whos % Define region to take EOF over eof_reg = [120 290 -60 60]; [xk, yk] = keep_var(eof_reg, lon, lat); ts = ts(:, yk, xk); lon = lon(xk); lat = lat(yk); [ntim, nlat, nlon] = size(ts); % Get rid of ice points oro = oro(:, yk, xk); oro = reshape(oro, ntim, nlat*nlon); oro = mean(oro); opt = find(oro == 0); %climmean = mean(clim); [ts, clim] = annave(ts); ts = reshape(ts, ntim, nlat, nlon); ts = cosweight(ts, lat); ts = reshape(ts, ntim, nlat*nlon); ts = ts(:, opt); ts = myrunning_ave(ts, 3); ts = myrunning_ave(ts, 3); c = ts*ts' / (nlat*nlon); [lam, pcs, per] = eof(c); pc10 = pcs(:,1:10); pc10 = (pc10 - ones(ntim, 1) * mean(pc10)) ./ (ones(ntim, 1) * std(pc10)); lds = pc10' * ts / ntim; lims = [0 360 -90 90] FRAME = eof_reg; num = 3; figure(num); sd(1); tem = NaN * ones(1, nlat*nlon); tem(1, opt) = lds(num, :); tem = reshape(tem, nlat, nlon); gcont(tem, [-1:.1:1]); dc xlabel(['Contour Interval = .1 K/std(pc', num2str(num), ')']); title(['T_S: EOF ' num2str(num) ', SOM']); sd(2) plot(pc10(:,num)); set(gca, 'XTick', 0:10:100, 'XTickLabel', 0:10:100); set(gca, 'YTick', -4:1:4); xlabel('Years'); ylabel('STD'); grid axis([0 100 -3 3]) cd /home/disk/tao/dvimont/matlab/CCM/SOM_run/SOM_Plots % Look at regressions against atmospheric data: cd /home/disk/hayes2/werner/matlab filin = '36yearly.cdf'; [precc, precl] = getnc(filin, 'PRECC', 'PRECL'); [lat, lon] = getll(filin); prec = (precc + precl) * 3600 * 24 * 1000; clear precc precl prec = getnc(filin, 'TS'); [lat, lon] = getll(filin); [ntim, nlat, nlon] = size(prec); prec = reshape(prec, ntim, nlat*nlon); [prec, clim] = annave(prec); clim = reshape(clim, 12, nlat, nlon); reg_pat = pc10' * prec / ntim; default_global num = 2; figure(num); sd(1); tem = NaN * ones(1, nlat*nlon); tem = reg_pat(num, :); tem = reshape(tem, nlat, nlon); gcont(tem, [-5:.1:5]); dc xlabel(['Contour Interval = .1 K/std(pc', num2str(num), ')']); title(['T_S: EOF ' num2str(num) ', SOM']); sd(2) plot(pc10(:,num)); set(gca, 'XTick', 0:10:100, 'XTickLabel', 0:10:100); set(gca, 'YTick', -4:1:4); xlabel('Years'); ylabel('STD'); grid axis([0 100 -3 3]) % Try this on only the Pacific: plims = [120 290 -30 70]; latk = find(lat >= plims(3) & lat <= plims(4)); lonk = find(lon >= plims(1) & lon <= plims(2)); tsnew = zeros(ntim, nlat*nlon); tsnew(:,opt) = ts; tsnew = reshape(tsnew, ntim, nlat, nlon); tsnew = tsnew(:, latk, lonk); [ptim, plat, plon] = size(tsnew); poro = reshape(oro(latk, lonk),1,plat*plon); popt = find(poro == 0); tsnew = reshape(tsnew, ntim, plat*plon); tsnew = tsnew(:, popt); c = tsnew'*tsnew ./ ntim; [plam, ploads, pper] = eof(c); plds = ploads(:,1:10) ./ (ones(length(popt),1) * std(ploads(:,1:10))); ppc = tsnew * plds ./ length(popt); plds = plds .* (ones(length(popt),1) * std(ppc)); ppc = ppc ./ (ones(ntim,1) * std(ppc)); num = 1; sd(1) tem = zeros(plat*plon, 1); tem(popt,1) = plds(:,num); tem = reshape(tem, plat, plon); gcont(-1*tem, plims, [-1:.1:1]); dc xlabel(['Contour Interval = .1 K/std(pc', num2str(num), ')']); title(['EOF ' num2str(num) ', SOM Pacific Domain']); sd(2) plot(-1*ppc(:,1)) set(gca, 'XTick', 120:120:1200, 'XTickLabel', 10:10:100); set(gca, 'YTick', -4:1:4); xlabel('Years'); grid % Get some atmospheric data that corresponds to some of this... nc = netcdf(filin, 'nowrite'); ps = nc{'PSL'}(:); nc = close(nc); ps = reshape(ps, ntim, nlat*nlon); [ps, pclim] = annave(ps); num = 1; pspat = pc' * ps ./ (100*ntim); global XAX YAX FRAME XAX = lon; YAX = lat; FRAME = lims; num = 1; sd(1) tem = reshape(pspat(num,:), nlat, nlon); %gcont(1*tem, lims, [-2:.25:2]); mcont(-1*tem, [-2:.25:2]); xlabel(['Contour Interval = 0.25 mb/std(pc', num2str(num), ')']); title(['SLP reg on PC ' num2str(num) ', SOM Global Domain']); %dc sd(2) plot(-1*pc(:,num)); set(gca, 'XTick', 120:120:1200, 'XTickLabel', 10:10:100); set(gca, 'YTick', -4:1:4); xlabel('Years'); ylabel('STD'); grid