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

load distinct_colors_20.mat


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


%----------------------------------
% Spatial composites over clusters 
%----------------------------------

% Load full SH data
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));

% define input and output variable
inputf='era5_sh_latlon.nc';
outputf=['era5_sh_pcp_cluster_comp.nc'];
% grab lat/lon from model, put it into the output file
eval(['!ncks -A -v latitude,longitude ',inputf,' ',outputf]);


%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)])

   % SAVE COMPOSITE
   cnm=['pcp_clust_' num2str(i) '_comp'];
   nccreate(outputf,cnm,'Dimensions',{'longitude',length(slon),'latitude',length(slat)});
   ncwrite(outputf,cnm,flip(cmp',2));

%   print(['/home/nsteiger/Dropbox/easterisle/ei_cluster_cmp' num2str(i) '_pcp.png'],'-dpng','-r300');
end

close all

% define input and output variable
inputf='era5_sh_latlon.nc';
outputf=['era5_sh_mslp_cluster_comp.nc'];
% grab lat/lon from model, put it into the output file
eval(['!ncks -A -v latitude,longitude ',inputf,' ',outputf]);

% 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)])

   % SAVE COMPOSITE
   cnm=['mslp_clust_' num2str(i) '_comp'];
   nccreate(outputf,cnm,'Dimensions',{'longitude',length(slon),'latitude',length(slat)});
   ncwrite(outputf,cnm,flip(cmp',2));

%   print(['/home/nsteiger/Dropbox/easterisle/ei_cluster_cmp' num2str(i) '_mslp.png'],'-dpng','-r300');
end

close all

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



%--------------------------------
% LAGGED COMPOSITES OF CLUSTER 4
%--------------------------------

% define input and output variable
inputf='era5_sh_latlon.nc';
outputf=['era5_sh_mslp_pcp_cluster4_comp_lags.nc'];
% grab lat/lon from model, put it into the output file
eval(['!ncks -A -v latitude,longitude ',inputf,' ',outputf]);

% MSLP
xf=reshape(xv_ds,length(slat),length(slon),length(dt));
xf2=reshape(xv2_ds,length(slat),length(slon),length(dt));

s=load('coast');mlat=s.lat;mlon=s.long;
lgs=[-4,-2,0,2,4];
%figure
for i=1:length(lgs)
   %subplot(length(lgs),1,i)
   figure
   hold on
   h=axesm('MapProjection','eqdcylin','MapLatLimit',[(min(slat)) (max(slat))],'MapLonLimit',[min(slon) max(slon)]);
   %cmp=mean(xf(:,:,cind),3)-xvmn;
   lix=rnc{4}+lgs(i);
   lix(lix<0)=[];% remove lags that go before beginning
   lix(lix>length(dt))=[];% remove lags that go after ending
   cmp=mean(xf(:,:,lix),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(lgs(i))])

   % SAVE COMPOSITE
   cnm=['mslp_clust_4_comp_lag_' num2str(lgs(i))];
   nccreate(outputf,cnm,'Dimensions',{'longitude',length(slon),'latitude',length(slat)});
   ncwrite(outputf,cnm,flip(cmp',2));

   cmp2=mean(xf2(:,:,lix),3);
   cnm=['pcp_clust_4_comp_lag_' num2str(lgs(i))];
   nccreate(outputf,cnm,'Dimensions',{'longitude',length(slon),'latitude',length(slat)});
   ncwrite(outputf,cnm,flip(cmp2',2));


%   print(['/home/nsteiger/Dropbox/easterisle/ei_cluster_cmp' num2str(i) '_mslp.png'],'-dpng','-r300');
end

close all




