Documentation of som_ts_eof

Global Index (short | long) | Local contents | Local Index (short | long)

Help text

filin = 'histape.cdf.36';

Cross-Reference Information

Listing of script som_ts_eof

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