
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
ot0=datetime(1955,1,1):datetime(2017,12,31);
ot=datenum(1955,1,1):datenum(2017,12,31);

% add NaNs for leap days since obs didn't include leap days...
leapdays=ot0(month(ot0)==2 & day(ot0)==29);
[~,ilp]=intersect(ot0,leapdays);
ida=setxor(1:length(ot0),ilp);
eirv=NaN(length(ot0),1);
eirv(ida)=eir(:);

% Remove leap days from time vector
%ot0(month(ot0)==2 & day(ot0)==29)=[];
%ot=datenum(ot0);

% ERA dates
dt=(datenum('01-Jan-1979'):datenum('31-Dec-2017'))';
dtt=(datetime('01-Jan-1979'):datetime('31-Dec-2017'))';

%pcp=ncread('/d1/nsteiger/climate-data/era5/era5_tp_swpac_daily_1979_2017.nc','tp');
pcp=ncread('/d1/nsteiger/climate-data/era5/era5_tp_sh_daily_1979_2017.nc','tp'); % larger area
pcp=flip(permute(pcp,[2 1 3]),1);
pcp(pcp<0)=0; % specify 0 for missing values

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

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

% 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,nlat,nlon]=nearest_latlon(slat,slon,eic_lat,eic_lon);


pcpv=reshape(pcp,length(slat)*length(slon),size(pcp,3));
mslpv=reshape(mslp,length(slat)*length(slon),size(mslp,3));

ei_p0=pcpv(nlm,:);
% In the obs, 0.1 mm is the lowest measurable amount, so impose that here for the analysis
ei_p0(ei_p0<1e-4)=0; % 1e-4 m = 0.1 mm
% get mslp at Easter Isl.
ei_m0=mslpv(nlm,:);

% monthly mean mslp
[ei_mslp_mon,mt2]=mean_ts(ei_m0,dt,'nanmean','month');

% Use obs for the rain, but replace missing obs with ERA data
[~,idte,idto]=intersect(dtt,ot0);
% where are missing values
ei_oe=eirv(idto); % obs + era5
ei_e0=ei_p0(idte)*1e3; % scale data to mm
ei_e0=round(ei_e0,1); % obs are rounded to nearest 0.1 mm, but model data is not...
nanobs=find(isnan(ei_oe));
ei_oe(nanobs)=ei_e0(nanobs); % replace nans with scaled era
ei_p=ei_oe; % use obs infilled with era5 % deseasoning makes negative rain!


% Seasonality of Easter Island precip
dseas=[1,2,10,11,12];
iseas=[3,8,9];
wseas=[4,5,6,7];


% mean of the annual precip
[mm_eir_p,ot2]=mean_ts(eirv,ot,'nansum','month'); % obs
mm_eir_p(mm_eir_p==0)=NaN;


% Reshape data
k=1;l=12;
mx_o=zeros(12,length(ot2)/12);
for n=1:length(ot2)/12
    mx_o(:,n)=mm_eir_p(k:l);
    k=k+12; l=l+12;
end

% Reshape mslp data
k=1;l=12;
mx_2=zeros(12,length(mt2)/12);
for n=1:length(mt2)/12
    mx_2(:,n)=ei_mslp_mon(k:l);
    k=k+12; l=l+12;
end



% How does seasonal data correlate with annual data?
[p_ei_yr,dt_yr]=mean_ts(eirv,ot,'nansum','year'); % obs
p_ei_yr(p_ei_yr==0)=NaN;

disp('corr(ann p, dseas)')
m1=mean(mx_o(dseas,:),1);
nancorr(m1,p_ei_yr)

disp('corr(ann p, tseas)')
m1=mean(mx_o(iseas,:),1);
nancorr(m1,p_ei_yr)

disp('corr(ann p, wseas)')
m1=mean(mx_o(wseas,:),1);
nancorr(m1,p_ei_yr)


% Average monthly precipitation
eio_mn=nanmean(mx_o,2);
eim_mn=nanmean(mx_2,2)./100; % MSLP

figure
plot(mean(eio_mn)*ones(12,1),'--','color',clrs.black,'linewidth',1.5)
hold on
plot(eio_mn,'-','color',clrs.grey,'linewidth',2)
% plot individual seasons by color
h1=plot(dseas,eio_mn(dseas),'s','color',clrs.orange,'linewidth',2,'markersize',8,'markerfacecolor',clrs.orange);
h2=plot(iseas,eio_mn(iseas),'o','color',clrs.black,'linewidth',2,'markersize',8,'markerfacecolor',clrs.black);
h3=plot(wseas,eio_mn(wseas),'^','color',clrs.green,'linewidth',2,'markersize',8,'markerfacecolor',clrs.green);
legend([h1,h2,h3],{'Dry season','Transition season','Wet season'})
ylabel('Average monthly 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)

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


figure
%plot(mean(eio_mn)*ones(12,1),'--','color',clrs.black,'linewidth',1.5)
hold on
plot(eim_mn,'-.','color',clrs.grey,'linewidth',2)
% plot individual seasons by color
h1=plot(dseas,eim_mn(dseas),'s','color',clrs.orange,'linewidth',2,'markersize',8,'markerfacecolor',clrs.orange);
h2=plot(iseas,eim_mn(iseas),'o','color',clrs.black,'linewidth',2,'markersize',8,'markerfacecolor',clrs.black);
h3=plot(wseas,eim_mn(wseas),'^','color',clrs.green,'linewidth',2,'markersize',8,'markerfacecolor',clrs.green);
legend([h1,h2,h3],{'Dry season','Transition season','Wet season'},'location','northwest')
ylabel('Average monthly MSLP (hPa)','fontsize',16)
box on
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)

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




% CALCULATE THE ANTICYCLONE FIELD ----------------------------------------------------------- 

% MEANS OF THE FIELD
% monthly mean mslp
[mslp_ds_m,mon_dt]=mean_ts(mslp,dt,'month');
% Annual mean 
[mn_mslp,tm_dt]=mean_ts(mslp,dt,'year');

% Reshape data
k=1;l=12; 
mX_o=zeros(size(mslp_ds_m,1),size(mslp_ds_m,2),12,length(mon_dt)/12);
for n=1:length(mon_dt)/12
    mX_o(:,:,:,n)=mslp_ds_m(:,:,k:l);
    k=k+12; l=l+12;
end


% % seasonal averages
% mslp_dry=squeeze(mean(mX_o(:,:,dseas,:),3));
% mslp_trn=squeeze(mean(mX_o(:,:,iseas,:),3));
% mslp_wet=squeeze(mean(mX_o(:,:,wseas,:),3));
 
% monthly climo
mslp_mclim=mean(mX_o,4);

% use edge of contour level ---------------------------------------
%pcut=1016; % pressure cut off (USE THIS FOR NCL PLOTTING PURPOSES)
%pcut=1019; % pressure cut off
pcut=1018; % pressure cut off
% reset lon values for this computation
lonmin=181;lonmax=288; % 221->220 deg, 288->287 deg = coastline of chile
latmin=41;latmax=76; % 41->50 degS, 76->15 degS

mslpX=mslp_mclim./100; %dd convert to mb
% find values greater than threshold to make a mask of regions
mslpX(mslpX<pcut)=0;
% just take region off South America
mslpX([1:latmin latmax:end],:,:)=0;
mslpX(:,[1:lonmin lonmax:end],:)=0;

% area-weight
W=cosd(repmat(slat,[1 length(slon)]));W=W(:);
mslpXV=reshape(mslpX,length(slat)*length(slon),12);
mslpXw=reshape(bsxfun(@times,mslpXV,W),length(slat),length(slon),12);


% center of mass location of these 1020 pressure blobs -------------------------
[Xs,Ys] = meshgrid(slon,slat);
cm=zeros(size(mslpX,3),2);
cmMass=zeros(size(mslpX,3),1);
for i=1:size(mslpX,3)
   [cm(i,:),cmMass(i)]=centroid(Xs(:),Ys(:),mslpXw(:,:,i));
end

% compute distance from easter island
cmD=dist_gc(eic_lat,eic_lon,cm(:,2),cm(:,1));

% Seasonal cycle of the anticyclone
figure
set(gcf,'defaultAxesColorOrder',[clrs.black; clrs.grey]);
yyaxis left
plot(1:12,360-cm(:,1),'--o','color',clrs.black,'linewidth',2)
ylabel('Pressure centroid longitude (deg W)','fontsize',16)
hold on
yyaxis right
plot(1:12,-1*cm(:,2),'-s','color',clrs.grey,'linewidth',2)
ylabel('Pressure centroid latitude (deg S)','fontsize',16)
xlim([1 12])
xticks(1:12)
xticklabels({'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'})
grid on
set(gca,'fontsize',14)

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


% normalized center of mass -----------------------

cmMn=cmMass./max(cmMass);

figure
hold on
plot(1:12,cmMn,'-o','color',clrs.black,'linewidth',2)
ylabel('Normalized pressure centroid mass','fontsize',16)
xlim([1 12])
xticks(1:12)
xticklabels({'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'})
grid on
box on
set(gca,'fontsize',14)

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


% location of centroid labeled by month -----------------------------------

% max mass
cmMx=max(mslpXV);
cmMxn=(rescale(cmMx)+0.5)*20; % rescaled in a format that's good for text font sizes

cmMmn=(rescale(cmMn)+0.5)*20; % rescaled in a format that's good for text font sizes

m={'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'};
elat=27.116667;
elon=109.366667;
figure
hold on
plot(360-cm(:,1),-1*cm(:,2),'.','color',clrs.white,'markersize',1)
%for i=1:12;text(360-cm(i,1),-1*cm(i,2),m{i},'fontsize',cmMmn(i),'color',clrs.black);end
for i=1:12;text(360-cm(i,1),-1*cm(i,2),m{i},'fontsize',cmMxn(i),'color',clrs.black);end
set(gca,'xdir','reverse','ydir','reverse')
xlabel('Pressure centroid longitude (deg W)','fontsize',16)
ylabel('Pressure centroid latitude (deg S)','fontsize',16)
plot(elon,elat,'p','color',[238 44 44]./255,'markersize',12)
box on
grid on
set(gca,'fontsize',14)

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

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

% distance of centroid away from easter island ---------------------------

figure
hold on
plot(1:12,cmD,'-','color',clrs.black,'linewidth',2)
for i=1:12;plot(i,cmD(i),'o','markersize',cmMxn(i),'color',clrs.black);end
ylabel('Distance from pressure centroid (km)','fontsize',16)
xlim([1 12])
xticks(1:12)
xticklabels({'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'})
grid on
box on
set(gca,'fontsize',14)

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



% % check contours
% figure
% for i=1:12
% pcolor(mslpX(:,:,i))
% print(['/home/nsteiger/Dropbox/easterisle/mslp_ann_' num2str(pcut) '_mon' num2str(i) '.png'],'-dpng','-r300');
% end

% save centroid vals as netcdf fields?
svncdf='n';
if svncdf=='y'

% define input and output variable
inputf='era5_sh_latlon.nc';
%outputf=['era5_mslp_monthlyclimo_1016contours_4nclplots.nc'];
outputf=['era5_mslp_monthlyclimo_1018contours_4centroidLocs.nc'];
%outputf=['era5_spac_1018contours_ann.nc'];

% grab lat/lon from model, put it into the output file
eval(['!ncks -A -v latitude,longitude ',inputf,' ',outputf]);

% Add time dimension
nccreate(outputf,'time','Dimensions',{'time',12});
ncwrite(outputf,'time',1:12);
ncwriteatt(outputf,'time','units','Monthly mean');
ncwriteatt(outputf,'time','long_name','Month');

% make april not have small blip of mslp near New Zealand
mslpX(:,[1:190 lonmax:end],4)=0;

for i=1:12

   Xv=mslpX(:,:,i);
   Xv=flip(Xv,1); % flip lats because python ncl screws it up in plotting for some reason...
   %cnm=['mslp1016_mon' num2str(i)];
   cnm=['mslp1018_mon' num2str(i)];
   %cnm=['mslp1018_' num2str(1978+i)];

   %nccreate(outputf,cnm,'Dimensions',{'longitude',length(slon),'latitude',length(slat)},'FillValue',1e20);
   nccreate(outputf,cnm,'Dimensions',{'longitude',length(slon),'latitude',length(slat)});
   ncwrite(outputf,cnm,permute(Xv,[2 1]));
%   ncwriteatt(outputf,cnm,'missing_value',1e20);

end
% add center of mass lat/lon values
nccreate(outputf,'cmlat','Dimensions',{'time',length(1:12)});
ncwrite(outputf,'cmlat',cm(:,2));
nccreate(outputf,'cmlon','Dimensions',{'time',length(1:12)});
ncwrite(outputf,'cmlon',cm(:,1));

end






return

% Plot observational seasonal cycle
figure
eirm=nanmean(eir,2);
plot(eirm,'color',clrs.grey,'linewidth',2)
hold on
plot(smooth(eirm,30,'lowess'),'r','linewidth',2)
ylabel('Daily mean precipitation (mm)','fontsize',16)
datetick('x','mmm')
grid on
%xticklabels({'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'})
set(gca,'fontsize',14)

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



% What season makes a dry vs wet year?
[mp,mdt]=mean_ts(ei_p,dt,'nansum','year');
[mp_obs,mdt_obs]=mean_ts(eirv,ot,'nansum','year');
[mp_mon,mdt_mon]=mean_ts(ei_p,dt,'nansum','month');

% % Reshape data
% k=1;l=12;
% mx_12=zeros(12,length(mdt_mon)/12);
% for n=1:length(mdt_mon)/12
%     mx_12(:,n)=mp_mon(k:l);
%     k=k+12; l=l+12;
% end

% Cluster analysis to find wet years
mp_obs(mp_obs==0)=NaN; % year of missing data...

[idx,cvs]=kmeans(mp_obs(:),3);
[~,i]=sort(cvs);
wids=find(idx==i(end)); % wet year indices
iids=find(idx==i(2)); % wet year indices
dids=find(idx==i(1)); % wet year indices

oids=cat(1,iids,dids); % just combine others as the 'null'
iwids=cat(1,iids,wids); % just combine others as the 'null'

%save('kmeans_wetdryint_yrs_easterisl.mat','wids','iids','dids','mp_obs','mdt_obs','idx','cvs')
%return


% % Which are wettest years?
% %[mps,mpi]=sort(mp,'descend');
% [mps,mpi]=sort(mp_obs,'descend');
% 
% 
% % Plot line color by dry/wet
% C=load('./colormaps/MPL_BrBG.mat');cmap_bg=C.cmap;
% %C=load('./colormaps/GMT_red2green.mat');%cmap_bg=C.cmap;
% %C=load('./colormaps/MPL_viridis.mat');
% % interpolate colormap onto small number of colors
% [X1,Y1] = meshgrid(1:3,1:length(C.cmap));
% [X,Y] = meshgrid(1:3,linspace(1,length(C.cmap),length(mpi)));
% cmapsm = interp2(X1,Y1,flip(C.cmap,1),X,Y);
% 
% 
% % looks like dry years are just missing wet events at any time of the year...
% %figure;for i=1:length(mpi); plot(1:12,mx_12(:,mpi(i)),'o','color',cmapsm(i,:),'linewidth',2);hold on;end
% figure;for i=1:length(mpi); plot(1:12,mx_o(:,mpi(i)),'o','color',cmapsm(i,:),'linewidth',2);hold on;end
% 
% % Top and bottom years
% %for i=[1:5 35:39]; plot(1:12,mx_12(:,mpi(i)),'-','color',cmapsm(i,:),'linewidth',2);hold on;end
% 
% % Overlay monthly mean
% %plot(mean(mx_12,2),'-o','color',clrs.black,'linewidth',2)
% %plot(mean(mx_12(:,mpi(1:13)),2),'-o','color',cmapsm(1,:),'linewidth',2)
% %plot(mean(mx_12(:,mpi(27:end)),2),'-o','color',cmapsm(end,:),'linewidth',2)
% 
% % 3rds + mean
% h1=plot(mean(mx_o(:,mpi(1:round(length(mpi)/4))),2),'--','color',cmapsm(1,:),'linewidth',2);
% h2=plot(nanmean(mx_o,2),'-','color',clrs.grey,'linewidth',2);
% h3=plot(mean(mx_o(:,mpi((end-1-round(length(mpi)/4)):end-1)),2),'-.','color',cmapsm(end,:),'linewidth',2);
% 
% % Overlay percentiles
% %mprc=prctile(mx_o',[25,50,75])';
% %plot(1:12,mprc(:,1),'-o','color',cmapsm(end,:),'linewidth',2)
% %plot(1:12,mprc(:,2),'-o','color',clrs.grey,'linewidth',2)
% %plot(1:12,mprc(:,3),'-o','color',cmapsm(1,:),'linewidth',2)
% 
% 
% legend([h1,h2,h3],{'Wettest 1/4','Average','Driest 1/4'})
% ylabel('Monthly precipitation (mm)','fontsize',16)
% grid on
% ax = gca;ax.XGrid='off';
% box on
% xlim([1 12])
% xticks(1:12)
% xticklabels({'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'})
% set(gca,'fontsize',14)
% cbh = colorbar ;
% cbh.Ticks = linspace(0,1,3);
% cbh.TickLabels = {'Dry','','Wet'};
% colormap(flip(cmapsm))
% %print(['/home/nsteiger/Dropbox/easterisle/ei_obs_wetdry_years_monthly_toppercent.png'],'-dpng','-r300');
% %print(['/home/nsteiger/Dropbox/easterisle/ei_obs_wetdry_years_monthly.png'],'-dpng','-r300');
% 


% Dry season, wet season, intermediate season

% 3rds + mean
wetyrs=mx_o(:,wids);
intyrs=mx_o(:,oids);
%dryyrs=mx_o(:,mpi((end+1-round(length(mpi)/3)):end));



% dry season
xd=[];gd=[];
x0=nansum(wetyrs(dseas,:),1);
%x0=nanmean(wetyrs(dseas,:),1);
g0=ones(length(x0),1);
xd=cat(1,xd(:),x0(:));
gd=cat(1,gd(:),g0(:));

x0=nansum(intyrs(dseas,:),1);
%x0=nanmean(intyrs(dseas,:),1);
g0=2*ones(length(x0),1);
xd=cat(1,xd(:),x0(:));
gd=cat(1,gd(:),g0(:));

% ks-test
%x0=nansum(wetyrs(dseas,:),1);
%x1=nansum(intyrs(dseas,:),1);
%[h,p]=kstest2(x0,x1)

% int season
x0=nansum(wetyrs(iseas,:),1);
%x0=nanmean(wetyrs(iseas,:),1);
g0=3*ones(length(x0),1);
xd=cat(1,xd(:),x0(:));
gd=cat(1,gd(:),g0(:));

x0=nansum(intyrs(iseas,:),1);
%x0=nanmean(intyrs(iseas,:),1);
g0=4*ones(length(x0),1);
xd=cat(1,xd(:),x0(:));
gd=cat(1,gd(:),g0(:));

% wet season

%xd(:,8)=nansum(dryyrs(wseas,:),1);
x0=nansum(wetyrs(wseas,:),1);
%x0=nanmean(wetyrs(wseas,:),1);
g0=5*ones(length(x0),1);
xd=cat(1,xd(:),x0(:));
gd=cat(1,gd(:),g0(:));

x0=nansum(intyrs(wseas,:),1);
%x0=nanmean(intyrs(wseas,:),1);
g0=6*ones(length(x0),1);
xd=cat(1,xd(:),x0(:));
gd=cat(1,gd(:),g0(:));

%xd(xd==0)=NaN;
%xd=xd(:);

%gd={1,2,3,1,2,3,1,2,3};
%gds=repmat(1:9,[21 ,1]);
%gd={1,2,1,2,1,2};
%gds=repmat(1:length(xd),[21 ,1]);


cg=[clrs.green;clrs.grey;clrs.green;clrs.grey;clrs.green;clrs.grey];%clrs.green];
%lbls={'Mid','Dry','Wet','Mid','Dry','Wet','Mid','Dry','Wet'};
%lbls={'Other','Wet','Other','Wet','Other','Wet'};
lbls={'Wet yrs','Other yrs','Wet yrs','Other yrs','Wet yrs','Other yrs'};

figure
%boxplot(xd,'colorgroup',gd,'symbol','','color',cg,'labels',lbls,'labelorientation','horizontal');
boxplot(xd,gd,'colorgroup',1:6,'symbol','','color',cg,'labels',lbls,'labelorientation','horizontal');
set(findobj(gca,'type','line'),'linew',2.5) % set linewidth thicker
grid on
ax = gca;ax.XGrid='off';
hold on
%scatter(gds(:),xd(:),'filled','MarkerFaceColor','k','MarkerFaceAlpha',0.4','jitter','on','jitterAmount',0.15)
scatter(gd(:),xd(:),'MarkerEdgeColor',clrs.grey,'MarkerEdgeAlpha',0.4','jitter','on','jitterAmount',0.15)
set(gca,'fontsize',14)
%xlabel('Year characteristic','fontsize',16)
ylabel('Accumulated precipitation (mm)','fontsize',16)

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

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



% Dry season, wet season, intermediate season

% 3rds + mean
dryyrs=mx_o(:,dids);
intyrs=mx_o(:,iwids);
%dryyrs=mx_o(:,mpi((end+1-round(length(mpi)/3)):end));



% dry season
xd=[];gd=[];
x0=nansum(dryyrs(dseas,:),1);
%x0=nanmean(wetyrs(dseas,:),1);
g0=ones(length(x0),1);
xd=cat(1,xd(:),x0(:));
gd=cat(1,gd(:),g0(:));

x0=nansum(intyrs(dseas,:),1);
%x0=nanmean(intyrs(dseas,:),1);
g0=2*ones(length(x0),1);
xd=cat(1,xd(:),x0(:));
gd=cat(1,gd(:),g0(:));

% ks-test
%x0=nansum(wetyrs(dseas,:),1);
%x1=nansum(intyrs(dseas,:),1);
%[h,p]=kstest2(x0,x1)

% int season
x0=nansum(dryyrs(iseas,:),1);
%x0=nanmean(wetyrs(iseas,:),1);
g0=3*ones(length(x0),1);
xd=cat(1,xd(:),x0(:));
gd=cat(1,gd(:),g0(:));

x0=nansum(intyrs(iseas,:),1);
%x0=nanmean(intyrs(iseas,:),1);
g0=4*ones(length(x0),1);
xd=cat(1,xd(:),x0(:));
gd=cat(1,gd(:),g0(:));

% wet season

%xd(:,8)=nansum(dryyrs(wseas,:),1);
x0=nansum(dryyrs(wseas,:),1);
%x0=nanmean(wetyrs(wseas,:),1);
g0=5*ones(length(x0),1);
xd=cat(1,xd(:),x0(:));
gd=cat(1,gd(:),g0(:));

x0=nansum(intyrs(wseas,:),1);
%x0=nanmean(intyrs(wseas,:),1);
g0=6*ones(length(x0),1);
xd=cat(1,xd(:),x0(:));
gd=cat(1,gd(:),g0(:));

%xd(xd==0)=NaN;
%xd=xd(:);

%gd={1,2,3,1,2,3,1,2,3};
%gds=repmat(1:9,[21 ,1]);
%gd={1,2,1,2,1,2};
%gds=repmat(1:length(xd),[21 ,1]);


cg=[clrs.brown;clrs.grey;clrs.brown;clrs.grey;clrs.brown;clrs.grey];%clrs.green];
%lbls={'Mid','Dry','Wet','Mid','Dry','Wet','Mid','Dry','Wet'};
%lbls={'Other','Wet','Other','Wet','Other','Wet'};
lbls={'Dry yrs','Other yrs','Dry yrs','Other yrs','Dry yrs','Other yrs'};

figure
%boxplot(xd,'colorgroup',gd,'symbol','','color',cg,'labels',lbls,'labelorientation','horizontal');
boxplot(xd,gd,'colorgroup',1:6,'symbol','','color',cg,'labels',lbls,'labelorientation','horizontal');
set(findobj(gca,'type','line'),'linew',2.5) % set linewidth thicker
grid on
ax = gca;ax.XGrid='off';
hold on
%scatter(gds(:),xd(:),'filled','MarkerFaceColor','k','MarkerFaceAlpha',0.4','jitter','on','jitterAmount',0.15)
scatter(gd(:),xd(:),'MarkerEdgeColor',clrs.grey,'MarkerEdgeAlpha',0.4','jitter','on','jitterAmount',0.15)
set(gca,'fontsize',14)
%xlabel('Year characteristic','fontsize',16)
ylabel('Accumulated precipitation (mm)','fontsize',16)

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

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








return

%--------------------
% Spatial composites 
%--------------------

clear slat slon t2m mslp pcp
% 
% lf_t2m='/d1/nsteiger/climate-data/era5/era5_t2m_scpac_daily_1979_2017.nc';
% %lf_t2m='/d1/nsteiger/climate-data/era5/era5_t2m_swpac_daily_1979_2017.nc';
lf_t2m='/d1/nsteiger/climate-data/era5/era5_t2m_sh_daily_1979_2017.nc';
% t2m=ncread(lf_t2m,'t2m'); % larger area
% t2m=flip(permute(t2m,[2 1 3]),1);
% 
slat=double(ncread(lf_t2m,'latitude'));
slat=flip(slat);
slon=double(ncread(lf_t2m,'longitude'));
% 
% %mslp=ncread('/d1/nsteiger/climate-data/era5/era5_msl_sh_daily_1979_2017.nc','msl'); % larger area
% %mslp=ncread('/d1/nsteiger/climate-data/era5/era5_msl_swpac_daily_1979_2017.nc','msl');
% mslp=ncread('/d1/nsteiger/climate-data/era5/era5_msl_scpac_daily_1979_2017.nc','msl'); % larger area
% mslp=flip(permute(mslp,[2 1 3]),1);
% 
pcp=ncread('/d1/nsteiger/climate-data/era5/era5_tp_sh_daily_1979_2017.nc','tp'); % larger area
% %pcp=ncread('/d1/nsteiger/climate-data/era5/era5_tp_swpac_daily_1979_2017.nc','tp');
%pcp=ncread('/d1/nsteiger/climate-data/era5/era5_tp_scpac_daily_1979_2017.nc','tp'); % larger area
pcp=flip(permute(pcp,[2 1 3]),1);
pcp(pcp<0)=0; % specify 0 for missing values

% Monthly mean of spatial fields
[pcp_mn,mdt_mon]=mean_ts(pcp,dt,'nansum','month');
pcp_amjj=mon2ann_sum(pcp_mn,4,7); % wet season
pcp_ondjf=mon2ann_sum(pcp_mn,10,2); % dry season


% CORRELATION MAPS
return

% Take annual average respecting water year...
pcp_ei_m2f=mon2ann_sum(mp_mon,3,2);
pcp_ei_amjj=mon2ann_sum(mp_mon,4,7);
pcp_ei_ondjf=mon2ann_sum(mp_mon,10,2);
% Sort annual values
%[~,is]=sort(pcp_ei_m2f,'descend');

xf=reshape(pcp_amjj,length(slat)*length(slon),size(pcp_amjj,3));
pcp_amjj_cr=zeros(length(slat)*length(slon),1);
pcp_amjj_pvl=zeros(length(slat)*length(slon),1);
for i=1:length(slat)*length(slon)
   %pcp_amjj_cr(i)=corr(pcp_ei_m2f(:),xf(i,1:end-1)');
   [pcp_amjj_cr(i),pcp_amjj_pvl(i)]=corr(pcp_ei_amjj(:),xf(i,:)');
end

xf=reshape(pcp_ondjf,length(slat)*length(slon),size(pcp_ondjf,3));
pcp_ondjf_cr=zeros(length(slat)*length(slon),1);
pcp_ondjf_pvl=zeros(length(slat)*length(slon),1);
for i=1:length(slat)*length(slon)
   %pcp_ondjf_cr(i)=corr(pcp_ei_m2f(:),xf(i,:)');
   [pcp_ondjf_cr(i),pcp_ondjf_pvl(i)]=corr(pcp_ei_ondjf(:),xf(i,:)');
end


% NaN-out high p-values
%ia=find(pcp_amjj_pvl>=0.05);
%io=find(pcp_ondjf_pvl>=0.05);
%pcp_amjj_cr(ia)=NaN;
%pcp_ondjf_cr(io)=NaN;




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

figure
s=load('coast');mlat=s.lat;mlon=s.long;
hold on
axesm('MapProjection','eqdcylin','MapLatLimit',[(min(slat)) (max(slat))],'MapLonLimit',[min(slon) max(slon)]);
h=pcolorm(slat,slon,reshape(pcp_amjj_cr,length(slat),length(slon)));
%set(h,'EdgeColor','none');
plotm(mlat,mlon,'color',[0.5 0.5 0.5])
caxis([-0.6 0.6])
colormap(cmap_bg)
colorbar
plotm(eic_lat,eic_lon,'p','color',clrs.red)
%set(gca,'xtick',[],'ytick',[])
title(['corr(E.I. AMJJ, AMJJ)'])
%print(['/home/nsteiger/Dropbox/easterisle/ei_corr_amjj_amjj.png'],'-dpng','-r300');


figure
s=load('coast');mlat=s.lat;mlon=s.long;
hold on
axesm('MapProjection','eqdcylin','MapLatLimit',[(min(slat)) (max(slat))],'MapLonLimit',[min(slon) max(slon)]);
h=pcolorm(slat,slon,reshape(pcp_ondjf_cr,length(slat),length(slon)));
%set(h,'EdgeColor','none');
plotm(mlat,mlon,'color',[0.5 0.5 0.5])
caxis([-0.6 0.6])
colormap(cmap_bg)
colorbar
plotm(eic_lat,eic_lon,'p','color',clrs.red)
%set(gca,'xtick',[],'ytick',[])
title(['corr(E.I. ONDJF, ONDJF)'])
%print(['/home/nsteiger/Dropbox/easterisle/ei_corr_ondjf_ondjf.png'],'-dpng','-r300');








return

%---------------------
% Compare against MJO
%---------------------
mjo0dt=(datetime('01-Jan-1979'):datetime('28-Aug-2018'))';
mjo=load('/home/nsteiger/Dropbox/data/OMI_MJO_index_NOAA.txt');
[~,~,midt]=intersect(dtt,mjo0dt); % use just up to 2017

mjoi=mjo(midt,7);
[mjoi_mon,~]=mean_ts(mjoi,midt,'nanmean','month');

pcpv=reshape(pcp,length(slat)*length(slon),size(pcp,3));
pcp_mjo_cr=zeros(length(slat)*length(slon),1);
for i=1:length(slat)*length(slon)
   pcp_mjo_cr(i)=corr(pcpv(i,:)',mjoi(:));
end



figure
s=load('coast');mlat=s.lat;mlon=s.long;
hold on
axesm('MapProjection','eqdcylin','MapLatLimit',[(min(slat)) (max(slat))],'MapLonLimit',[min(slon) max(slon)]);
h=pcolorm(slat,slon,reshape(pcp_mjo_cr,length(slat),length(slon)));
%set(h,'EdgeColor','none');
plotm(mlat,mlon,'color',[0.5 0.5 0.5])
%caxis([-0.6 0.6])
colormap(cmap_bg)
colorbar
plotm(eic_lat,eic_lon,'p','color',clrs.red)
%set(gca,'xtick',[],'ytick',[])
title(['corr(MJO, PCP)'])
%print(['/home/nsteiger/Dropbox/easterisle/ei_corr_ondjf_ondjf.png'],'-dpng','-r300');


