%enso and station data
loaddap('http://iridl.ldeo.columbia.edu/expert/SOURCES/.Indices/.nino/.EXTENDED/.NINO34/T/3/runningAverage/T/%28Sep-Nov%29VALUES/T/%281960%29%282009%29RANGE/dods')

loaddap('http://iridl.ldeo.columbia.edu/expert/SOURCES/.Indices/.nino/.EXTENDED/.NINO34/T/3/runningAverage/T/%28Dec-Feb%29VALUES/T/%281960%29%282010%29RANGE/dods')


loaddap('http://iridl.ldeo.columbia.edu/expert/ds:/SOURCES/.NOAA/.NCEP/.CPC/.CAMS/.station/.precipitation/IWMO/91348/91353/91366/91371/91376/91378/91425/91502/91503/91507/91520/91533/91601/91610/91623/91629/91631/91636/91643/91724/91800/91925/VALUES/a:/.lat/:a:/.lon/:a/:ds/dods')


loaddap('http://iridl.ldeo.columbia.edu/expert/SOURCES/.NOAA/.NCEP/.CPC/.CAMS/.station/.precipitation/IWMO/91348/91353/91366/91371/91376/91378/91425/91502/91503/91507/91520/91533/91601/91610/91623/91629/91631/91636/91643/91724/91800/91925/VALUES/.prcp/T/3/runningAverage/T/%28Dec-Feb%29VALUES/T/%28Dec%201960%29%28Dec%202008%29RANGE/dods')

loaddap('http://iridl.ldeo.columbia.edu/expert/SOURCES/.NOAA/.NCEP/.CPC/.CAMS/.station/.precipitation/IWMO/97014/VALUES/.prcp/T/3/runningAverage/T/%28Dec-Feb%29VALUES/T/%28Dec%201960%29%28Dec%202008%29RANGE/dods')

obs = prcp.prcp;
corr2(obs(~isnan(obs)),v(~isnan(obs),1)')  
'Menado, Indonesia'

iok = find(mean(~isnan(obs),2)==1);

lat=lat.lat;
lon=lon.lon;
[lat(iok(1)),lon(iok(1))]
[lat(iok(2)),lon(iok(2))]

ttl1='Majuro Island in the Marshall Islands';
ttl2='Atuona, French Polynesia';

loaddap('http://iridl.ldeo.columbia.edu/expert/home/.tippett/.scratch/.HadSST2/.ssta/T/3/runningAverage/T/%28Sep-Nov%29VALUES/T/%281960%29%282009%29RANGE/lat/-20/20/RANGE/lon/130/260/RANGE/dods')

X = lon;
Y = lat;

ssta = ssta.ssta;
msk = ~isnan(mean(ssta,3));

[ny,nx,nt] = size(ssta);

x = reshape(ssta,ny*nx,nt);
xm = mean(x,2);
xa = x - xm(:,ones(1,nt)); 
xa(~msk,:)=0;

[u,s,v] = svd(xa,'econ');
u(~msk,:)=NaN;
my_pcolor(X,Y,reshape(u(:,1),ny,nx))

y=obs';
b = regress(y,[ones(nt,1),v(:,1)]);
yhat = [ones(nt,1),v(:,1)]*b;
sqrt(nanmean((yhat-y).^2))

b_glm = glmfit(v(:,1),y,'poisson');
yhat_glm = glmval(b_glm,v(:,1),'log');
sqrt(nanmean((yhat_glm-y).^2))

xx = [1.25*min(v(:,1)):0.01:1.2*max(v(:,1))]';
yline = [ones(size(xx)),xx]*b;

ycurve = glmval(b_glm,xx,'log');

plot(v(:,1)',y,'.')
hold on
plot(xx,[yline,ycurve]);
hold off



