
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 daily precipitation (mm)','fontsize',16)
% datetick('x','mmm')
% set(gca,'fontsize',14)
% 
% %print(['/home/nsteiger/Dropbox/easterisle/ei_rain_avg.png'],'-dpng','-r300');



disp('Load ERA5 data')


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

%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_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

% % load 0.125 deg resolution precip data from ERA5
% pcp_hr=ncread('/home/nsteiger/Dropbox/data/era5_tp_easterisle_daily_1979_1981.nc','tp');
% pcp_hr=flip(permute(pcp_hr,[2 1 3]),1);
% pcp_hr(pcp_hr<0)=0; % specify 0 for missing values
% 
% hlat=double(ncread('/home/nsteiger/Dropbox/data/era5_tp_easterisle_daily_1979_1981.nc','latitude'));
% hlat=flip(hlat);
% hlon=double(ncread('/home/nsteiger/Dropbox/data/era5_tp_easterisle_daily_1979_1981.nc','longitude'));

%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=double(ncread('/d1/nsteiger/climate-data/era5/era5_t2m_scpac_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_scpac_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,nlat,nlon]=nearest_latlon(slat,slon,eic_lat,eic_lon);
%[nlm_h,nlat_h,nlon_h]=nearest_latlon(hlat,hlon,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;

%pcpvh=reshape(pcp_hr,length(hlat)*length(hlon),size(pcp_hr,3));
pcpv=reshape(pcp,length(slat)*length(slon),size(pcp,3));
mslpv=reshape(mslp,length(slat)*length(slon),size(mslp,3));
t2mv=reshape(t2m,length(slat)*length(slon),size(t2m,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
ei_m0=mslpv(nlm,:);
ei_t0=t2mv(nlm,:);


%eih_p0=pcpvh(nlm_h,:);
% In the obs, 0.1 mm is the lowest measurable amount, so impose that here for the analysis
%eih_p0(eih_p0<1e-4)=0; % 1e-4 m = 0.1 mm





% PLOT THE OBSERVATIONS AGAINST ERA5

% 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'))';
dsht=(datetime('01-Jan-1979'):datetime('31-Dec-1981'))';

% mean of the annual precip
[mm_eir_p,ot2]=mean_ts(eirv,ot,'nansum','month'); % obs
mm_eir_p(mm_eir_p==0)=NaN;
[mm_eie_p,et2]=mean_ts(ei_p0(:)*1e3,dt,'nansum','month'); % era5



% 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 data
k=1;l=12;
mx_e=zeros(12,length(et2)/12);
for n=1:length(et2)/12
    mx_e(:,n)=mm_eie_p(k:l);
    k=k+12; l=l+12;
end

% Average monthly precipitation
eio_mn=nanmean(mx_o,2);
eie_mn=nanmean(mx_e,2);


% figure
% %plot(eie_mn,'-o','color',clrs.purple,'linewidth',2)
% %hold on
% plot(eio_mn,'-o','color',clrs.purple,'linewidth',2)
% ylabel('Mean 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');
% 
% % 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('Mean daily 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');




% Plot obs and ERA5 on overlaping interval
[~,~,idt0]=intersect(dsht,dtt);
[~,~,idt]=intersect(dsht,ot0);

% figure
% plot(eirv(idt),'o-','color',clrs.blue)
% hold on
% plot(eih_p0*1e3,'o-','color',clrs.red)
% plot(ei_p0(idt0)*1e3,'o-','color',clrs.purple)
% legend('Obs','ERA5_{hr}','ERA5')


% 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


% Remove seasonal cycle of MSLP
%ei_p=ei_p0; % deseasoning makes negative rain!
ei_p=ei_oe; % use obs infilled with era5 % deseasoning makes negative rain!
%ei_p=deseason(dt,ei_p0);
%ei_m=deseason(dt,ei_m0);
%ei_t=deseason(dt,ei_t0);

% Combine data for multivariate clustering
% NOTE: DOING THIS SEEMS TO MUDDY UP THE CLUSTERS BECAUSE MSLP TENDS TO DOMINATE
%Xd(:,1)=rescale(ei_p);
%Xd(:,2)=rescale(ei_m);
%Xd(:,3)=rescale(ei_t);

% Normalize data?
%Xd=normalize(Xd,'range');

%=======
% ERA5
%=======

%x=1:length(ei_p(:));
%eipv0=eip(:);eipv=eipv0;
% Pull out no-rain
%eip(eip<1e-4)=0; % 1e-4 m = 0.1 mm
%nrn=find(eipv0<1e-4);
%eipv(nrn)=NaN;



%-------------------
% K-MEANS CLUSTERING
%-------------------

% Multivariate clustering (doesn't seem to work as well, gets drug around by MSLP)
%[cidx,cent]=kmeans(Xd,4,'distance','cityblock');

% Just cluster on precip
%[cidx,cent]=kmeans(ei_p(:),4,'distance','cityblock','MaxIter',1000);
% Load pre-computed k-means for consistency (sometimes algorithm struggles to converge)
load('kmeans_clusters_8-May-2019.mat')


% Sort by amount of precip
[cent_s,is]=sort(cent);


%rnc{1}=nrn; % no rain
rnc{1}=find(cidx==is(1)); 
rnc{2}=find(cidx==is(2)); 
rnc{3}=find(cidx==is(3)); 
rnc{4}=find(cidx==is(4)); 
%rnc4=find(cidx==4);


% Plot representative years of data
exdt=(datetime('01-Jan-1979'):datetime('31-Dec-1981'))';
[~,~,exid]=intersect(exdt,dtt);


figure
%plot(dtt,ei_p,'.') 
hold on
plot(dtt(rnc{1}),ei_p(rnc{1}),'.','color',clrs.black)
plot(dtt(rnc{2}),ei_p(rnc{2}),'o','color',clrs.red)
plot(dtt(rnc{3}),ei_p(rnc{3}),'^','color',clrs.blue)
plot(dtt(rnc{4}),ei_p(rnc{4}),'s','color',clrs.green)

xlim([dtt(exid(1)) dtt(exid(end))])
%xlim([dtt(1) dtt(end)])
datetick('x','yyyy')
box on
set(gca,'fontsize',16)
ylabel('Daily precipitation (mm)','fontsize',18)
xlabel('Year','fontsize',18)
legend('Cluster 1','Cluster 2','Cluster 3','Cluster 4','location','northoutside','orientation','horizontal')

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

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


%============================
% Box plot of clusters
%============================


lbls={'1','2','3','4'};
%grp={1,2,3,4,1,2,3,4};
%grp={1,2,3,4};
grp = [zeros(length(rnc{1}),1); ones(length(rnc{2}),1); 2*ones(length(rnc{3}),1); 3*ones(length(rnc{4}),1)];
cgrp=[clrs.black;clrs.red;clrs.blue;clrs.green];
cdt=[ei_p(rnc{1});ei_p(rnc{2});ei_p(rnc{3});ei_p(rnc{4})];

figure
%boxplot(varSplt,'colorgroup',grp,'PlotStyle','compact','symbol','','color',cgrp,...
%  'labels',lbls,'labelorientation','horizontal');

boxplot(cdt,grp,'colorgroup',grp,'symbol','o','color',cgrp,...
  'labels',lbls)%,'labelorientation','horizontal');
%boxplot(cdt,grp,'colorgroup',grp,'symbol','','color',cgrp,...
%  'labels',lbls)%,'labelorientation','horizontal');
set(findobj(gca,'type','line'),'linew',2) % set linewidth thicker
%title(seas)
%ylabel('Precipitation','fontsize',16,'fontweight','bold')
ylabel('Daily precipitation (mm)','fontsize',18)
xlabel('Cluster','fontsize',18)
%ylim([-0.5 0.95])
grid on
ax = gca;ax.XGrid='off';
%set(gca,'fontsize',22,'fontweight','bold')
set(gca,'fontsize',16)
%set(findobj(gca,'Type','text'),'FontSize',18,'VerticalAlignment','top') % labels done differently in boxplot

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

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






% mean annual precipitation
[mp,mdt]=mean_ts(ei_p,dt,'nansum','year');

% get the mean of each cluster by replacing non-cluster values by nans and nansumming
pc1=NaN(size(ei_p));
pc2=NaN(size(ei_p));
pc3=NaN(size(ei_p));
pc4=NaN(size(ei_p));

pc1(rnc{1})=ei_p(rnc{1});
pc2(rnc{2})=ei_p(rnc{2});
pc3(rnc{3})=ei_p(rnc{3});
pc4(rnc{4})=ei_p(rnc{4});

[mpc1,mdt]=mean_ts(pc1,dt,'nansum','year');
[mpc2,mdt]=mean_ts(pc2,dt,'nansum','year');
[mpc3,mdt]=mean_ts(pc3,dt,'nansum','year');
[mpc4,mdt]=mean_ts(pc4,dt,'nansum','year');

% NaN-out entirely ERA-based year
[Y,~,~,~,~,~]=datevec(mdt);
ni=find(Y==1990);
mpc1(ni)=NaN;
mpc2(ni)=NaN;
mpc3(ni)=NaN;
mpc4(ni)=NaN;
mp(ni)=NaN;

% Total amount of rain
figure
hold on
plot(Y,mp,'-*','color',clrs.purple,'linewidth',2)
plot(Y,mpc1,'-.','color',clrs.black,'linewidth',2)
plot(Y,mpc2,'-o','color',clrs.red,'linewidth',2)
plot(Y,mpc3,'-^','color',clrs.blue,'linewidth',2)
plot(Y,mpc4,'-s','color',clrs.green,'linewidth',2)
legend('Total','Cluster 1','Cluster 2','Cluster 3','Cluster 4','location','northoutside','orientation','horizontal')
%datetick('x','yyyy')
ylabel('Annual precipitation (mm)','fontsize',16)
xlim([Y(1) Y(end)])
xlabel('Year','fontsize',16)
%xtickangle(45)
grid on
box on
%ax = gca;ax.XGrid='off';
set(gca,'fontsize',14)

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

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

disp(['Cluster 1 avg = ' num2str(nanmean(mpc1))])
disp(['Cluster 2 avg = ' num2str(nanmean(mpc2))])
disp(['Cluster 3 avg = ' num2str(nanmean(mpc3))])
disp(['Cluster 4 avg = ' num2str(nanmean(mpc4))])
disp(['Total avg = ' num2str(nanmean(mp))])


disp(['corr(total,C1) = ' num2str(nancorr(mp,mpc1))])
disp(['corr(total,C2) = ' num2str(nancorr(mp,mpc2))])
disp(['corr(total,C3) = ' num2str(nancorr(mp,mpc3))])
disp(['corr(total,C4) = ' num2str(nancorr(mp,mpc4))])



% % What season makes a dry vs wet year?
% [mp,mdt]=mean_ts(ei_p,dt,'nansum','year');
% [mp_mon,mdt]=mean_ts(ei_p,dt,'nansum','month');
% 
% % Reshape data
% k=1;l=12;
% mx_12=zeros(12,length(et2)/12);
% for n=1:length(et2)/12
%     mx_12(:,n)=mp_mon(k:l);
%     k=k+12; l=l+12;
% end
% mxnm=bsxfun(@minus,mx_12,mean(mx_12,1));
% 
% % Which are wettest years?
% [~,mpi]=sort(mp,'descend');
% 
% 
% % Plot line color by dry/wet
% C=load('./colormaps/MPL_BrBG.mat');cmap_bg=C.cmap;
% % 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
% 
% % 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
% 
% 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 = {'Driest','','Wettest'};
% 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');
% 




%mp_ann=mean(mean_ts(ei_p,dt,'nansum','year'));
%fr_ei_p=ei_p/ma_ei_p; % fraction of average mean annual precipitation


% Count up the frequencies of each cluster per year
[y0,m0,~,~,~,~] = datevec(dt);
CompArray = datenum(y0,1,1);
ut = unique(CompArray);
fpc=zeros(4,length(ut));
for k = 1:length(ut)
   fpc(1,k)=sum(~isnan(pc1(CompArray==ut(k))));
   fpc(2,k)=sum(~isnan(pc2(CompArray==ut(k))));
   fpc(3,k)=sum(~isnan(pc3(CompArray==ut(k))));
   fpc(4,k)=sum(~isnan(pc4(CompArray==ut(k))));
end

fpc(:,ni)=NaN;

% Total amount of rain
figure
hold on
plot(Y,fpc(1,:),'-.','color',clrs.black,'linewidth',2)
plot(Y,fpc(2,:),'-o','color',clrs.red,'linewidth',2)
plot(Y,fpc(3,:),'-^','color',clrs.blue,'linewidth',2)
plot(Y,fpc(4,:),'-s','color',clrs.green,'linewidth',2)
legend('Cluster 1','Cluster 2','Cluster 3','Cluster 4','location','northoutside','orientation','horizontal')
%datetick('x','yyyy')
ylabel('Precipitation cluster frequency','fontsize',16)
xlim([Y(1) Y(end)])
xlabel('Year','fontsize',16)
%xtickangle(45)
grid on
box on
%ax = gca;ax.XGrid='off';
set(gca,'fontsize',14)

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

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

[avgpc4,mdt]=mean_ts(pc4,dt,'nanmean','year');
disp('-------------------------------------')
disp('Avg c4 precip vs freq of c4 precip')
disp(['corr(total,Frq C4) = ' num2str(nancorr(mp(:),fpc(4,:)))])
disp(['corr(total,Avg C4) = ' num2str(nancorr(mp(:),avgpc4(:)))])
disp(['std(Avg C4) = ' num2str(std(avgpc4))])
disp('-------------------------------------')

disp(['corr(total,Frq C1) = ' num2str(nancorr(mp(:),fpc(1,:)))])
disp(['corr(total,Frq C2) = ' num2str(nancorr(mp(:),fpc(2,:)))])
disp(['corr(total,Frq C3) = ' num2str(nancorr(mp(:),fpc(3,:)))])
disp(['corr(total,Frq C4) = ' num2str(nancorr(mp(:),fpc(4,:)))])

disp(['corr(pc1,Frq C1) = ' num2str(nancorr(mpc1(:),fpc(1,:)))])
disp(['corr(pc2,Frq C2) = ' num2str(nancorr(mpc2(:),fpc(2,:)))])
disp(['corr(pc3,Frq C3) = ' num2str(nancorr(mpc3(:),fpc(3,:)))])
disp(['corr(pc4,Frq C4) = ' num2str(nancorr(mpc4(:),fpc(4,:)))])


return


% Count up the average monthly frequencies
CompArray = datenum(y0,m0,1);
ut = unique(CompArray);
fmpc=zeros(4,length(ut));
for k = 1:length(ut)
   fmpc(1,k)=sum(~isnan(pc1(CompArray==ut(k))));
   fmpc(2,k)=sum(~isnan(pc2(CompArray==ut(k))));
   fmpc(3,k)=sum(~isnan(pc3(CompArray==ut(k))));
   fmpc(4,k)=sum(~isnan(pc4(CompArray==ut(k))));
end

fmpcr=reshape(fmpc,4,12,length(ut)/12);
mf=mean(fmpcr,3);

figure
hold on
plot(mf(1,:),'-.','color',clrs.black,'linewidth',2)
plot(mf(2,:),'-o','color',clrs.red,'linewidth',2)
plot(mf(3,:),'-^','color',clrs.blue,'linewidth',2)
plot(mf(4,:),'-s','color',clrs.green,'linewidth',2)

legend('Cluster 1','Cluster 2','Cluster 3','Cluster 4','location','northoutside','orientation','horizontal')
ylabel('Mean frequency','fontsize',16)
grid on
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)

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


% BOX PLOT OF MEAN FREQUENCIES

lbls={'1','2','3','4'};
%grp={1,2,3,4,1,2,3,4};
%grp={1,2,3,4};
grp = [zeros(length(mdt),1); ones(length(mdt),1); 2*ones(length(mdt),1); 3*ones(length(mdt),1)];
cgrp=[clrs.black;clrs.red;clrs.blue;clrs.green];
cdt=fpc';

figure
%boxplot(varSplt,'colorgroup',grp,'PlotStyle','compact','symbol','','color',cgrp,...
%  'labels',lbls,'labelorientation','horizontal');

boxplot(cdt(:),grp,'colorgroup',grp,'symbol','o','color',cgrp,...
  'labels',lbls)%,'labelorientation','horizontal');
%boxplot(cdt,grp,'colorgroup',grp,'symbol','','color',cgrp,...
%  'labels',lbls)%,'labelorientation','horizontal');
set(findobj(gca,'type','line'),'linew',2) % set linewidth thicker
%title(seas)
%ylabel('Precipitation','fontsize',16,'fontweight','bold')
ylabel('Annual frequency','fontsize',18)
xlabel('Cluster','fontsize',18)
%ylim([-0.5 0.95])
grid on
ax = gca;ax.XGrid='off';
%set(gca,'fontsize',22,'fontweight','bold')
set(gca,'fontsize',16)
%set(findobj(gca,'Type','text'),'FontSize',18,'VerticalAlignment','top') % labels done differently in boxplot

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

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


mfnm=bsxfun(@minus,mf,mean(mf,2));

disp('Average events per month for each cluster')
mean(mf,2)

figure
hold on
plot(mfnm(1,:),'-.','color',clrs.black,'linewidth',2)
plot(mfnm(2,:),'-o','color',clrs.red,'linewidth',2)
plot(mfnm(3,:),'-^','color',clrs.blue,'linewidth',2)
plot(mfnm(4,:),'-s','color',clrs.green,'linewidth',2)

legend('Cluster 1','Cluster 2','Cluster 3','Cluster 4','location','northoutside','orientation','horizontal')
ylabel('Mean frequency anomaly','fontsize',18)
grid on
box on
xlim([1 12])
xticks(1:12)
xticklabels({'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'})
set(gca,'fontsize',16)
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0.5 2 7.25 4.5]; % [left bottom width height]

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


% Mean of clusters by month
[pc1_mon,mondt]=mean_ts(pc1,dt,'nansum','month');
[pc2_mon,mondt]=mean_ts(pc2,dt,'nansum','month');
[pc3_mon,mondt]=mean_ts(pc3,dt,'nansum','month');
[pc4_mon,mondt]=mean_ts(pc4,dt,'nansum','month');


% Reshape data
k=1;l=12;
pc1m_12=zeros(12,length(mondt)/12);
pc2m_12=zeros(12,length(mondt)/12);
pc3m_12=zeros(12,length(mondt)/12);
pc4m_12=zeros(12,length(mondt)/12);
for n=1:length(mondt)/12
    pc1m_12(:,n)=pc1_mon(k:l);
    pc2m_12(:,n)=pc2_mon(k:l);
    pc3m_12(:,n)=pc3_mon(k:l);
    pc4m_12(:,n)=pc4_mon(k:l);
    k=k+12; l=l+12;
end


figure
hold on
plot(mean(pc1m_12,2),'-.','color',clrs.black,'linewidth',2)
plot(mean(pc2m_12,2),'-o','color',clrs.red,'linewidth',2)
plot(mean(pc3m_12,2),'-^','color',clrs.blue,'linewidth',2)
plot(mean(pc4m_12,2),'-s','color',clrs.green,'linewidth',2)

legend('Cluster 1','Cluster 2','Cluster 3','Cluster 4','location','northoutside','orientation','horizontal')
ylabel('Mean precipitation (mm)','fontsize',18)
grid on
box on
xlim([1 12])
xticks(1:12)
xticklabels({'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'})
set(gca,'fontsize',16)
fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0.5 2 7.25 4.5]; % [left bottom width height]

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


return


% %----------------------------------
% % Spatial composites over clusters 
% %----------------------------------
% 
% % Load full SH data
% 
% clear slat slon t2m mslp pcp
% 
% t2m=ncread('/d1/nsteiger/climate-data/era5/era5_t2m_sh_daily_1979_2017.nc','t2m'); % larger area
% t2m=flip(permute(t2m,[2 1 3]),1);
% 
% 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);
% 
% 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
% 
% 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_sh_daily_1979_2017.nc','longitude'));
% 
% 
% 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));
% xv3=reshape(t2m,length(slat)*length(slon),length(dt));
% xv_ds=zeros(size(xv));
% xv2_ds=zeros(size(xv));
% xv3_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,:));
%    xv3_ds(i,:)=deseason(dt,xv3(i,:));
% end
% 
% 
% 
% 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);
% 
% 
% % 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
% s=load('coast');mlat=s.lat;mlon=s.long;
% for i=1:length(rnc)
%    figure
%    %subplot(length(rnc),1,i)
%    hold on
%    h=axesm('MapProjection','eqdcylin','MapLatLimit',[(min(slat)) (max(slat))],'MapLonLimit',[min(slon) max(slon)]);
%    %cmp=mean(xf(:,:,cind),3)-xvmn;
%    cmp=mean(xf(:,:,rnc{i}),3);
%    %h=contour(X,Y,cmp);
%    h=pcolorm(slat,slon,cmp);
%    %set(h,'EdgeColor','none');
%    plotm(mlat,mlon,'color',[0.5 0.5 0.5])
%    caxis([-5e-3 5e-3]) % Make the axis a log-axis?
%    colormap(cmap_bg)
%    plotm(eic_lat,eic_lon,'p','color',clrs.red)
%    %set(gca,'xtick',[],'ytick',[])
%    title(['Precip. ' num2str(i)])
% 
% %   print(['/home/nsteiger/Dropbox/easterisle/ei_cluster_cmp' num2str(i) '_pcp.png'],'-dpng','-r300');
% end
% 
% 
% 
% % MSLP
% xf=reshape(xv_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
% s=load('coast');mlat=s.lat;mlon=s.long;
% for i=1:length(rnc)
%    figure
%    %subplot(length(rnc),1,i)
%    hold on
%    h=axesm('MapProjection','eqdcylin','MapLatLimit',[(min(slat)) (max(slat))],'MapLonLimit',[min(slon) max(slon)]);
%    %cmp=mean(xf(:,:,cind),3)-xvmn;
%    cmp=mean(xf(:,:,rnc{i}),3);
%    %h=contour(X,Y,cmp);
%    h=pcolorm(slat,slon,cmp);
%    %set(h,'EdgeColor','none');
%    plotm(mlat,mlon,'color',[0.5 0.5 0.5])
%    caxis([-3e2 3e2]) % Make the axis a log-axis?
%    colormap(cmap_br)
%    plotm(eic_lat,eic_lon,'p','color',clrs.red)
%    %set(gca,'xtick',[],'ytick',[])
%    title(['MSLP ' num2str(i)])
% 
% %   print(['/home/nsteiger/Dropbox/easterisle/ei_cluster_cmp' num2str(i) '_mslp.png'],'-dpng','-r300');
% end
% 
% 
% % T2M
% xf=reshape(xv3_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
% s=load('coast');mlat=s.lat;mlon=s.long;
% for i=1:length(rnc)
%    figure
%    %subplot(length(rnc),1,i)
%    hold on
%    h=axesm('MapProjection','eqdcylin','MapLatLimit',[(min(slat)) (max(slat))],'MapLonLimit',[min(slon) max(slon)]);
%    %cmp=mean(xf(:,:,cind),3)-xvmn;
%    cmp=mean(xf(:,:,rnc{i}),3);
%    %h=contour(X,Y,cmp);
%    h=pcolorm(slat,slon,cmp);
%    %set(h,'EdgeColor','none');
%    plotm(mlat,mlon,'color',[0.5 0.5 0.5])
%    caxis([-0.6 0.6]) % Make the axis a log-axis?
%    colormap(cmap_br)
%    plotm(eic_lat,eic_lon,'p','color',clrs.red)
%    %set(gca,'xtick',[],'ytick',[])
%    title(['T2m ' num2str(i)])
% 
% %   print(['/home/nsteiger/Dropbox/easterisle/ei_cluster_cmp' num2str(i) '_t2m.png'],'-dpng','-r300');
% end


% %------------------------------------------
% % Plot some of the underlying spatial data
% 
% 
% %xf=reshape(xv_ds,length(slat),length(slon),length(dt));
% xf=reshape(xv2_ds,length(slat),length(slon),length(dt));
% 
% figure
% for i=1:25
%    subplot(5,5,i)
% %   h=contour(X,Y,xf(:,:,rnc{4}(i)));
% %   caxis([-3e2 3e2]) % Make the axis a log-axis?
% %   colormap(cmap_br)
%    h=pcolor(X,Y,xf(:,:,rnc{4}(i)));
%    set(h, 'EdgeColor', 'none');
%    caxis([0 0.2]) % Make the axis a log-axis?
%    %colormap(cmap_bg)
%    hold on
%    plot(eic_lon,eic_lat,'kp')
%    set(gca,'xtick',[],'ytick',[])
%    title(num2str(i))
% end
% suptitle('MSLP')
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 9 6.5]; % [left bottom width height]











%==============
% SOM analysis
%==============

%load('output_som30_era5_mslp_26-Mar-2019_20:04:25.mat','bmus')
%load('output_som35_era5_mslp_26-Mar-2019_20:08:26.mat','bmus')
%load('output_som40_era5_mslp_26-Mar-2019_21:24:30.mat','bmus')

%load('output_som40_era5_mslp_23-May-2019_16:42:35.mat','bmus','somRows','somCols')
%load('output_som48_era5_mslp_23-May-2019_16:26:27.mat','bmus','somRows','somCols')
%load('output_som54_era5_mslp_23-May-2019_16:42:18.mat','bmus','somRows','somCols')

% standardized data, smaller area, more soms
load('output_som54_era5_mslp_27-May-2019_21:42:16.mat','bmus','somRows','somCols');
%load('output_som48_era5_mslp_27-May-2019_21:40:32.mat','bmus','somRows','somCols');


% SOM frequency vs SOMs during clusters

figure
histogram(bmus,somRows*somCols,'normalization','pdf','facecolor',clrs.grey) 
hold on
histogram(bmus(rnc{1}),somRows*somCols,'normalization','pdf')
title('Cluster 1')
ylim([0 0.07])
ylabel('Frequency (norm)')
xlabel('SOM')
set(gca,'fontsize',14)
%print(['/home/nsteiger/Dropbox/easterisle/ei_cluster_hist1.png'],'-dpng','-r300');

figure
histogram(bmus,somRows*somCols,'normalization','pdf','facecolor',clrs.grey) 
hold on
histogram(bmus(rnc{2}),somRows*somCols,'normalization','pdf')
title('Cluster 2')
ylim([0 0.07])
ylabel('Frequency (norm)')
xlabel('SOM')
set(gca,'fontsize',14)
%print(['/home/nsteiger/Dropbox/easterisle/ei_cluster_hist2.png'],'-dpng','-r300');

figure
histogram(bmus,somRows*somCols,'normalization','pdf','facecolor',clrs.grey) 
hold on
histogram(bmus(rnc{3}),somRows*somCols,'normalization','pdf')
title('Cluster 3')
ylim([0 0.07])
ylabel('Frequency (norm)')
xlabel('SOM')
set(gca,'fontsize',14)
%print(['/home/nsteiger/Dropbox/easterisle/ei_cluster_hist3.png'],'-dpng','-r300');


figure
histogram(bmus,somRows*somCols,'normalization','pdf','facecolor',clrs.grey) 
hold on
histogram(bmus(rnc{4}),somRows*somCols,'normalization','pdf')
title('Cluster 4')
ylim([0 0.07])
ylabel('Frequency (norm)')
xlabel('SOM')
set(gca,'fontsize',14)
%print(['/home/nsteiger/Dropbox/easterisle/ei_cluster_hist4.png'],'-dpng','-r300');






% %----------------------
% % Compare against ENSO
% 
% 
% [X_obs,xlat,xlon,time0] = load_obs_Smon('bearth');
% X_obs=reshape(X_obs,length(xlat),length(xlon),length(time0)*12);
% [t_nino0]=nino(X_obs,xlat,xlon,3,'n');
% [~,ier,ibe] = intersect(1979:2017,time0);
% tninoann=mon2ann(t_nino0,1,12);
% nancorr(tninoann(ibe),mpc4(ier))
% 
% 
% % East pacific ITCZ
% [mpcp,mdt]=mean_ts(pcp,dt,'nansum','year');
% [itcz_mx,~] = itcz_lon(mpcp(41:end,:,:),slat(41:end),slon,[130 260]);
% nancorr(itcz_mx,mpc4)
% 
% 
% %------------------------------
% % Compare against SAM
% %------------------------------
% 
% % Load SH MSLP
% mslp_sh=ncread('/d1/nsteiger/climate-data/era5/era5_msl_sh_daily_1979_2017.nc','msl'); % larger area
% mslp_sh=flip(permute(mslp_sh,[2 1 3]),1);
% 
% shlat=double(ncread('/d1/nsteiger/climate-data/era5/era5_t2m_sh_daily_1979_2017.nc','latitude'));
% shlat=flip(shlat);
% shlon=double(ncread('/d1/nsteiger/climate-data/era5/era5_t2m_sh_daily_1979_2017.nc','longitude'));
% 
% % extract mslp at 45 and 60 degrees south
% lt1=-65;idl_65=find((abs(shlat-lt1)-min(abs(shlat-lt1)))<1e-5);
% lt1=-40;idl_40=find((abs(shlat-lt1)-min(abs(shlat-lt1)))<1e-5);
% 
% % monthly mean mslp
% [mslp_ds_m,mon_dt]=mean_ts(reshape(mslp_sh,length(shlat),length(shlon),length(dt)),dt,'month');
% 
% zm_mslp_40=mean(squeeze(mslp_ds_m(idl_40,:,:)),1);
% zm_mslp_65=mean(squeeze(mslp_ds_m(idl_65,:,:)),1);
% 
% sami=zscore(zm_mslp_40)-zscore(zm_mslp_65);
% 
% [pc4_mon,~]=mean_ts(pc4,dt,'nansum','month');
% 
% 
% [pc_t,~]=mean_ts(ei_p,dt,'nansum','month');
% 
% 
% % correlation map of SAM index with MSLP field
% 
% [pcp_mon,~]=mean_ts(pcp,dt,'nansum','month');
% pcpv_mon=reshape(pcp_mon,length(slat)*length(slon),length(mon_dt));
% mslp_dsmv=reshape(mslp_ds_m,length(slat)*length(slon),length(mon_dt));
% 
% cr_sam_mslp=zeros(length(slat)*length(slon),1);
% cr_sam_pcp=zeros(length(slat)*length(slon),1);
% for i=1:length(slat)*length(slon)
%    cr_sam_mslp(i)=corr(mslp_dsmv(i,:)',sami(:));
%    cr_sam_pcp(i)=corr(pcpv_mon(i,:)',sami(:));
% end
% 
% figure
% pcolor(reshape(cr_sam_mslp,length(slat),length(slon)))
% title('MSLP')
% 
% figure
% pcolor(reshape(cr_sam_pcp,length(slat),length(slon)))
% title('PCP')
% 





 %--------------------------------------------------
 % Composite over the big rain vs small rain events
 %--------------------------------------------------
% 
% % Remove seasonal cycle of MSLP
% disp('Removing the seasonal cycle...')
% dt=(datenum('01-Jan-1979'):datenum('31-Dec-2017'))';
% xv=reshape(mslp,length(slat)*length(slon),length(dt));
% xv_ds=zeros(size(xv));
% parfor i=1:length(slat)*length(slon)
%    xv_ds(i,:)=deseason(dt,xv(i,:));
% end
% 
% % Just take ocean values
% xvf=NaN(length(slat)*length(slon),length(dt));
% xvf(oids,:)=xv_ds(oids,:);
% xf=reshape(xvf,length(slat),length(slon),length(dt));
% 
% 
% C=load('./colormaps/NCV_blue_red.mat');cmap_br=C.cmap;
% C=load('./colormaps/MPL_BrBG.mat');cmap_bg=C.cmap;
% s=load('coast');mlat=s.lat;mlon=s.long;
% 
% 
% [X,Y] = meshgrid(slon,slat);
% 
% for i=1:length(rnc)
% 
%    figure
%    %axesm('mapprojection','eqdcylin')
%    %h=axesm('MapProjection','wagner4','MapLatLimit',[(min(slat)) (max(slat))],'MapLonLimit',[min(slon) max(slon)]);
%    hold on
%    cmp=mean(xf(:,:,rnc{i}),3);
%    h=contourf(X,Y,cmp,-5e2:20:5e2);
%    %h=contourfm(slat,slon,cmp);
%    %set(h,'EdgeColor','none');
%    %caxis([-5e-3 5e-3]) % Make the axis a log-axis?
%    caxis([-5e2 5e2]) % Make the axis a log-axis?
%    colormap(cmap_br)
%    colorbar
%    plot(eic_lon,eic_lat,'gp')
%    %set(gca,'xtick',[],'ytick',[])
%    title(['MSLP cluster ' num2str(i)])
% 
%    fig = gcf;fig.PaperUnits = 'inches';
%    fig.PaperPosition = [0 2 7 5]; % [left bottom width height]
%    %print(['/home/nsteiger/Dropbox/easterisle/ei_cmp_mslp_cluster' num2str(i) '.png'],'-dpng','-r300');
% 
% end






% %==============================
% % COMPARISON WITH OBSERVATIONS
% %==============================
% 
% 
% % Matlab's quartile outlier detector does a good job separating the two scales of rain
% tf=isoutlier(eip(:),'quartiles');
% bgrn=find(tf==1); % big rain
% norn=find(eip(:)==0); % no rain
% smrn=setdiff(1:length(eip),union(bgrn,norn)); % small rain
% 
% 
% % Plot a sample
% x=1:length(eip(:));
% figure;
% plot(x,eip(:),'.') 
% xlim([x(1) x(720)])
% hold on
% plot(x(bgrn),eip(bgrn),'sr')
% plot(x(norn),eip(norn),'om')
% xlabel('Day')
% ylabel('Precip')
% legend('Rain','Outlier','No rain')
% title('ERA5')
% %print(['/home/nsteiger/Dropbox/easterisle/ei_rain_outlier.png'],'-dpng','-r300');
% 
% 
% % OBSERVATIONS
% % Matlab's quartile outlier detector does a good job separating the two scales of rain
% tf_o=isoutlier(eir(:),'quartiles');
% bgrn_o=find(tf_o==1); % big rain
% norn_o=find(eir(:)==0); % no rain
% smrn_o=setdiff(1:length(eir(:)),union(bgrn_o,norn_o)); % small rain
% 
% % Plot a sample
% x=1:length(eir(:));
% figure;
% plot(x,eir(:),'.') 
% xlim([x(1) x(720)])
% hold on
% plot(x(bgrn_o),eir(bgrn_o),'sr')
% plot(x(norn_o),eir(norn_o),'om')
% xlabel('Day')
% ylabel('Precip')
% legend('Rain','Outlier','No rain')
% title('Observations')
% %print(['/home/nsteiger/Dropbox/easterisle/ei_obs_rain_outlier.png'],'-dpng','-r300');
% 
% 
% 
% 
% %-------------------
% % Cluster analysis
% %-------------------
% eirv0=eir(:);eirv=eirv0;
% % Pull out no-rain
% nrn=find(eirv0<0.1);
% eirv(nrn)=NaN;
% 
% 
% %cidx=kmeans(eirv(:),4);
% %cidx=kmeans(eirv(:),4,'distance','cityblock');
% cidx=kmeans(eirv(:),3,'distance','cityblock');
% %cidx=kmeans(eirv(:),2,'distance','cityblock');
% 
% rnc1=find(cidx==1); 
% rnc2=find(cidx==2); 
% rnc3=find(cidx==3); 
% %rnc4=find(cidx==4);
% rnc4=nrn; % no rain
% 
% figure
% plot(x,eir(:),'.') 
% xlim([x(1) x(720)])
% hold on
% plot(x(rnc1),eirv0(rnc1),'sr')
% plot(x(rnc2),eirv0(rnc2),'om')
% plot(x(rnc3),eirv0(rnc3),'ok')
% plot(x(rnc4),eirv0(rnc4),'pg')
% xlabel('Day')
% ylabel('Precip')
% legend('Rain','C1','C2','C3','No rain')
% %legend('Rain','C1','C2','No rain')
% %legend('Rain','C1','C2','C3')
% title('Observations')
% %print(['/home/nsteiger/Dropbox/easterisle/ei_obs_rain_cluster3_norain.png'],'-dpng','-r300');
% 















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











