
tic
clear

load distinct_colors_20.mat

% Load observatoinal data
eir0 = xlsread('/home/nsteiger/Dropbox/data/easter_island_rain.xlsx');
eir_yr=eir0(1,:);
eir=eir0(2:end,:);
% Start and end dates of data
%tt=datetime(1955,1,1):datetime(2017,12,31);
% Remove leap days
%tt(month(tt)==2 & day(tt)==29) = [];
%eirt=timetable(tt,eir(:)); % doesn't seem to work...

figure
%mdt=month(datetime(1955,1,1):datetime(1955,12,31));
eirm=nanmean(eir,2);
plot(eirm,'color',clrs.grey,'linewidth',2)
hold on
plot(smooth(eirm,30,'lowess'),'r','linewidth',2)
ylabel('Mean precipitation (mm)','fontsize',16)
datetick('x','mmm')
set(gca,'fontsize',14)

%print(['/home/nsteiger/Dropbox/easterisle/ei_rain_avg.png'],'-dpng','-r300');


% Plot distributions?




disp('Load ECHAM5-WISO data')


t2m=ncread('/d3/nsteiger/echam5wiso/t2m_mon_1871_2011.nc','temp2');
t2m=flip(permute(t2m,[2 1 3]),1);

%slp=ncread('/d3/nsteiger/echam5wiso/slp_mon_1871_2011.nc','slp');
%slp=flip(permute(slp,[2 1 3]),1);

pcp=ncread('/d3/nsteiger/echam5wiso/precip_mon_1871_2011.nc','precip');
pcp=flip(permute(pcp,[2 1 3]),1);

%wnd=ncread('/d3/nsteiger/echam5wiso/wind10_mon_1871_2011.nc','wind10');
%wnd=flip(permute(wnd,[2 1 3]),1);

dDp=squeeze(ncread('/d3/nsteiger/echam5wiso/dD_precip_mon_1871_2011.nc','wisoaprt_d'));
dDp=flip(permute(dDp,[2 1 3]),1);


lat=ncread('/d3/nsteiger/echam5wiso/slp_mon_1871_2011.nc','lat');
lat=flip(lat);
lon=ncread('/d3/nsteiger/echam5wiso/slp_mon_1871_2011.nc','lon');
nlns=find(lon<0);lon(nlns)=lon(nlns)+360;


% Easter island lats/lons: -27.116667, -109.366667
eic_lat=-27.116667;
eic_lon=360-109.366667;

% Where in the model is Easter Island?
nlm=nearest_latlon(lat,lon,eic_lat,eic_lon);

% What area to analyze? 
% lat(25)=-62.24, lat(70)=-11.77, lon(185)=207, lon(252)=282 
%et=30:72;en=185:252;
et=25:103;en=125:259; % broad tropical Pacific box
eilat=lat(et);eilon=lon(en);

% Where is ocean?
load('/d3/nsteiger/echam5wiso/sftlf_fx_ECHAM5-wiso.mat')
omsk=sftlf;omsk(isnan(omsk))=1;omsk(omsk==100)=0;


% Does ECHAM have realistic precipitation?
v0=reshape(pcp,length(lat)*length(lon),size(t2m,3));
eie_pcp_mon=v0(nlm,:);
% convert precip rate to mm/day
eie_mn=mean(reshape(86400*eie_pcp_mon,12,141),2);

% cludgy way of taking a monthly mean...
ei_rmon(1)=mean(eirm(1:31));ei_rmon(2)=mean(eirm(32:59));
ei_rmon(3)=mean(eirm(60:90));ei_rmon(4)=mean(eirm(91:120));
ei_rmon(5)=mean(eirm(121:151));ei_rmon(6)=mean(eirm(152:181));
ei_rmon(7)=mean(eirm(182:212));ei_rmon(8)=mean(eirm(213:243));
ei_rmon(9)=mean(eirm(244:273));ei_rmon(10)=mean(eirm(274:304));
ei_rmon(11)=mean(eirm(305:334));ei_rmon(12)=mean(eirm(335:365));

figure
plot(eie_mn,'-o','color',clrs.purple,'linewidth',2)
hold on
plot(ei_rmon,'-o','color',clrs.red,'linewidth',2)
ylabel('Mean precipitation (mm)','fontsize',16)
grid on
xlim([1 12])
xticks(1:12)
xticklabels({'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'})
set(gca,'fontsize',14)
legend('ECHAM5-wiso','Obs')
%print(['/home/nsteiger/Dropbox/easterisle/ei_rain_avg_vs_echam.png'],'-dpng','-r300');


% Sesonal averaging
mon_avg_o=4;mon_avg_f=3;
%mon_avg_o=6;mon_avg_f=8;
%mon_avg_o=12;mon_avg_f=2;

%Xvmon=slp(et,en,:);
%Xvmon=pcp(et,en,:);
Xvmon=t2m(et,en,:);
pomsk=omsk(et,en);
oids=find(pomsk(:)==1);
W=cosd(repmat(eilat,[1 length(eilon)]));W=W(:);
wt=W(oids);

Xv=mon2ann(Xvmon,mon_avg_o,mon_avg_f);
xv=reshape(Xv,length(et)*length(en),size(Xv,3));

% Just take ocean values
xv=xv(oids,:);

%------------------
%  Calculate SOMs
%------------------

% Detrend (not really much of a trend at all in the data...)
D=detrend(xv')';
%D=xv;

% Area weighting
%W=cosd(repmat(eilat,[1 length(eilon)]));W=W(:);


toc

somRows=2;somCols=2;
disp(['SOMs = ' num2str(somRows*somCols)])
init='l'; % (l)inear or (r)andom initialization
% D[lat*lon,time]
disp('Computing SOM...')
[ptrns0,bmus,rmsErr,pfrq,sMap,sD] = som_make_njs(D,somRows,somCols,init,wt);
disp('SOM completed');toc

% Expand out the ptrns
ptrns=NaN(somRows*somCols,length(et)*length(en));
ptrns(:,oids)=ptrns0;


%disp('Saving output...')

%timestp2=strrep(datestr(now),' ','_');

%save(['output_som' num2str(somRows*somCols) '_yrs' num2str(analyrs0(1)) num2str(analyrs0(end)) '_' timestp2 '.mat'],'-v7.3')

toc

%----------------
% PLOT MSLP SOMS
%----------------

C=load('./colormaps/NCV_blue_red.mat');cmap_br=C.cmap;
C=load('./colormaps/MPL_BrBG.mat');cmap_bg=C.cmap;

[X,Y] = meshgrid(eilon,eilat);


%----------------
% PLOT MSLP SOMS
%----------------
% figure
% for i=1:somRows*somCols
%    subplot(somRows,somCols,i)
%    h=contour(X,Y,reshape(ptrns0(i,:),length(et),length(en)));
%    caxis([-1 1])
%    colormap(cmap_br)
%    hold on
%    plot(eic_lon,eic_lat,'rp')
% end
% suptitle('MSLP')

%----------------
% PLOT T2M SOMS
%----------------

figure
for i=1:somRows*somCols
   subplot(somRows,somCols,i)
   %h=contour(X,Y,reshape(ptrns0(i,:),length(et),length(en)));
   h=pcolor(X,Y,reshape(ptrns(i,:),length(et),length(en)));
   set(h,'EdgeColor','none');
   caxis([-1.5 1.5])
   colormap(cmap_br)
   hold on
   plot(eic_lon,eic_lat,'kp')
end
suptitle('T2m')
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 9 5]; % [left bottom width height]

%print(['/home/nsteiger/Dropbox/easterisle/ei_soms' num2str(somRows*somCols) '_sst.png'],'-dpng','-r300');

%---------------------------
% COMPOSITE OTHER VARIABLES
%---------------------------

% Precip
Xvmon=pcp(et,en,:);
%Xvmon=slp(et,en,:);
Xv=mon2ann(Xvmon,mon_avg_o,mon_avg_f);
xvmn=mean(Xv,3);

figure
for i=1:somRows*somCols
   subplot(somRows,somCols,i)
   cind=find(bmus==i);
   cmp=mean(Xv(:,:,cind),3)-xvmn;
   %h=contour(X,Y,cmp);
   h=pcolor(X,Y,cmp);
   set(h,'EdgeColor','none');
   caxis([-7e-6 7e-6]) % Make the axis a log-axis?
   colormap(cmap_bg)
   hold on
   plot(eic_lon,eic_lat,'rp')
end
suptitle('\Delta PCP')

fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 9 5]; % [left bottom width height]

%print(['/home/nsteiger/Dropbox/easterisle/ei_cmp' num2str(somRows*somCols) '_pcp.png'],'-dpng','-r300');

% dD
Xvmon=dDp(et,en,:);
%Xvmon=slp(et,en,:);
Xv=mon2ann(Xvmon,mon_avg_o,mon_avg_f);
xvmn=mean(Xv,3);

figure
for i=1:somRows*somCols
   subplot(somRows,somCols,i)
   cind=find(bmus==i);
   cmp=mean(Xv(:,:,cind),3)-xvmn;
   %h=contour(X,Y,cmp);
   h=pcolor(X,Y,cmp);
   set(h,'EdgeColor','none');
   caxis([-5 5]) 
   colormap(cmap_br)
   hold on
   plot(eic_lon,eic_lat,'kp')
end
suptitle('\delta D precip')

fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 9 5]; % [left bottom width height]

%print(['/home/nsteiger/Dropbox/easterisle/ei_cmp' num2str(somRows*somCols) '_dDp.png'],'-dpng','-r300');


% What does PDSI do?
load('/d3/nsteiger/echam5wiso/pdsi_output_AWC_c_7_u2_gcm_14-Feb-2019_19:29:29.mat','pdsi_lndonly')
%e_pdsi=reshape(pdsi_lndonly,160,320,1692);
e_pdsi=pdsi_lndonly;

% Grap PDSI over EI
ei_pdsi=mon2ann(e_pdsi(nlm,:),mon_avg_o,mon_avg_f);

v0=reshape(t2m,length(lat)*length(lon),size(t2m,3));
ei_t2m=mon2ann(v0(nlm,:),mon_avg_o,mon_avg_f);

v0=reshape(dDp,length(lat)*length(lon),size(t2m,3));
ei_dDp=mon2ann(v0(nlm,:),mon_avg_o,mon_avg_f);

v0=reshape(pcp,length(lat)*length(lon),size(t2m,3));
ei_pcp_mon=v0(nlm,:);
ei_pcp=mon2ann(v0(nlm,:),mon_avg_o,mon_avg_f);


% How are different variables related?
figure
scatter(ei_dDp,ei_pdsi)

figure;scatter(ei_dDp,ei_pdsi)

figure;scatter(ei_dDp,ei_t2m)

figure;scatter(ei_dDp,ei_pcp)
















