
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...
%
% % Plot observational seasonal cycle
% 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 ERA5 data')


t2m=ncread('/d1/nsteiger/climate-data/era5/era5_t2m_swpac_daily_1979_2017.nc','t2m');
t2m=flip(permute(t2m,[2 1 3]),1);

mslp=ncread('/d1/nsteiger/climate-data/era5/era5_msl_swpac_daily_1979_2017.nc','msl');
mslp=flip(permute(mslp,[2 1 3]),1);

pcp=ncread('/d1/nsteiger/climate-data/era5/era5_tp_swpac_daily_1979_2017.nc','tp');
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);


slat=double(ncread('/d1/nsteiger/climate-data/era5/era5_t2m_swpac_daily_1979_2017.nc','latitude'));
slat=flip(slat);
slon=double(ncread('/d1/nsteiger/climate-data/era5/era5_t2m_swpac_daily_1979_2017.nc','longitude'));

% Where is ocean?
lat=ncread('/d1/nsteiger/climate-data/era5/landseamask_orography_era5_1deg.nc','latitude');
lat=flip(lat);
lon=ncread('/d1/nsteiger/climate-data/era5/landseamask_orography_era5_1deg.nc','longitude');
[~,islt]=intersect(lat,slat);
[~,isln]=intersect(lon,slon);

lsm=ncread('/d1/nsteiger/climate-data/era5/landseamask_orography_era5_1deg.nc','lsm');
lsm=flip(permute(lsm,[2 1]),1);
lmsk=lsm(islt,isln);
ln=find(lmsk>=0.5);oc=find(lmsk<0.5);
omsk=lmsk;omsk(ln)=0;omsk(oc)=1;

% Weight the data over the ocean by area
oids=find(omsk(:)==1);
W=cosd(repmat(slat,[1 length(slon)]));W=W(:);
wt=W(oids);

% 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(slat,slon,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;eilon=lon;



% % 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');


disp('Removing the seasonal cycle...')
dt=(datenum('01-Jan-1979'):datenum('31-Dec-2017'))';
xv=reshape(mslp,length(slat)*length(slon),length(dt));
xv2=reshape(pcp,length(slat)*length(slon),length(dt));
%xv=reshape(t2m,length(slat)*length(slon),length(dt));
xv_ds=zeros(size(xv));
xv2_ds=zeros(size(xv));
parfor i=1:length(slat)*length(slon)
   xv_ds(i,:)=deseason(dt,xv(i,:));
   xv2_ds(i,:)=deseason(dt,xv2(i,:));
end

% Just take ocean values
xvo=xv_ds(oids,:);
%xvo=xv(oids,:);disp('Data is not deseasoned...')


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

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

% Just take a subset for faster processing...
%D=D(:,1:10*366);
%xv2_ds=xv2_ds(:,1:10*366);

toc

somRows=6;somCols=5;
disp(['SOMs = ' num2str(somRows*somCols)])
init='r'; % (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(slat)*length(slon));
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')

%----------------
% 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(slon,slat);


%----------------
% PLOT MSLP SOMS
%----------------
figure
for i=1:somRows*somCols
   subplot(somRows,somCols,i)
   h=contour(X,Y,reshape(ptrns(i,:),length(slat),length(slon)));
   caxis([-1.5 1.5])
   colormap(cmap_br)
   hold on
   plot(eic_lon,eic_lat,'kp')
   set(gca,'xtick',[],'ytick',[])
end
suptitle('MSLP')
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 9 5]; % [left bottom width height]

timestp2=strrep(datestr(now),' ','_');
print(['/home/nsteiger/Dropbox/easterisle/ei_soms' num2str(somRows*somCols) '_mslp_' timestp2 '.png'],'-dpng','-r300');

%----------------
% 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
xf=reshape(xv2_ds,length(slat),length(slon),length(dt));
%xf=reshape(xv2_ds,length(slat),length(slon),length(1:10*366));
xvmn=squeeze(mean(xf,3));

figure
for i=1:somRows*somCols
   subplot(somRows,somCols,i)
   cind=find(bmus==i);
   %cmp=mean(xf(:,:,cind),3)-xvmn;
   cmp=mean(xf(:,:,cind),3);
   %h=contour(X,Y,cmp);
   h=pcolor(X,Y,cmp);
   set(h,'EdgeColor','none');
   caxis([-5e-3 5e-3]) % Make the axis a log-axis?
   colormap(cmap_bg)
   hold on
   plot(eic_lon,eic_lat,'kp')
   set(gca,'xtick',[],'ytick',[])
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_' timestp2 '.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)
















