
clear


disp('Load ERA5 SOM data')


% RAW DATA
%load('output_som12_era5_mslp_nostd_02-Apr-2019_13:56:42.mat')
%load('output_som16_era5_mslp_nostd_02-Apr-2019_13:51:54.mat')
%load('output_som20_era5_mslp_nostd_02-Apr-2019_13:48:19.mat')
%load('output_som25_era5_mslp_nostd_02-Apr-2019_12:18:30.mat')
%load('output_som30_era5_mslp_nostd_02-Apr-2019_12:12:33.mat')
%load('output_som35_era5_mslp_nostd_02-Apr-2019_12:24:14.mat')

% STANDARDIZED DATA
%load('output_som25_era5_mslp_26-Mar-2019_21:18:26.mat')
%load('output_som30_era5_mslp_26-Mar-2019_20:04:25.mat') 
%load('output_som35_era5_mslp_26-Mar-2019_20:08:26.mat')
%load('output_som40_era5_mslp_26-Mar-2019_21:24:30.mat')

% On clustered data
%load('output_som12_era5_mslp_nostd_13-May-2019_20:00:11.mat')

% standardized data, larger geographical area
%load('output_som40_era5_mslp_23-May-2019_16:42:35.mat');
%load('output_som48_era5_mslp_23-May-2019_16:26:27.mat');
%load('output_som54_era5_mslp_23-May-2019_16:42:18.mat');

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

%----------------
% 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])
%   caxis([-2000 2000])
   colormap(cmap_br)
   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]

% %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([-4e-3 4e-3]) % 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('\Delta PCP')

fig = gcf;fig.PaperUnits = 'inches';
fig.PaperPosition = [0 2 9 6.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)



% % eip => Easter Island Precip in ERA5
% % eir => Easter Island Precip in obs
% % In the obs, 0.1 mm is the lowest measurable amount, so impose that here for the analysis
% eip(eip<1e-4)=0; % 1e-4 m = 0.1 mm
% % Make zero rain NaNs for the clustering algorithm
% eipv=eip(:);
% nrn=find(eipv<1e-4);
% eipv(nrn)=NaN;
% 
% 
% %-------------------
% % K-MEANS CLUSTERING
% %-------------------
% 
% % Assume the number of clusters
% [cidx,cent]=kmeans(eipv(:),3,'distance','cityblock');
% 
% % Sort by amount of precip
% [cent_s,is]=sort(cent);
% 
% % Which days belong to which cluster?
% rnc{1}=nrn; % no rain
% rnc{2}=find(cidx==is(1)); 
% rnc{3}=find(cidx==is(2)); 
% rnc{4}=find(cidx==is(3)); 
% 
% 
% % Deseasoned data:
% % xv_ds => mslp, xv2_ds => pcp
% 
% % Just take ocean values
% %xvo=xv_ds(oids,:);
% 
% 
% % Understand state transitions
% tM=transM(bmus,1); % transition matrix
% mc = dtmc(tM); % make a markov chain object
% 
% % PLOT A DIRECTED GRAPH
% %figure;graphplot(mc) % plot the full directed graph 





% % Plot states corresponding to SOM nodes
% C=load('./colormaps/NCV_blue_red.mat');cmap_br=C.cmap;
% C=load('./colormaps/MPL_BrBG.mat');cmap_bg=C.cmap;
% 
% % What does the normalized data that goes into the SOM look like?
% addpath('/home/nsteiger/MATLAB/somtoolbox_2005')
% Dn=som_normalize(D','var')';
% % Expand out the data fields
% Df=NaN(size(xv_ds));
% Df(oids,:)=Dn;
% 
% % Let's look at SOM node 30
% b30=find(bmus==30);
% 
% [X,Y] = meshgrid(slon,slat);
% figure
% for i=1:9
%    subplot(3,3,i)
%    h=contour(X,Y,reshape(Df(:,b30(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',[])
%    title(num2str(i))
% end
% suptitle('Standardized MSLP')
% 
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 8 6]; % [left bottom width height]
% %print(['/home/nsteiger/Dropbox/easterisle/ei_underlyingdata_soms' num2str(somRows*somCols) '_mslp_std_' timestp2 '.png'],'-dpng','-r300');
% 
% figure
% for i=1:9
%    subplot(3,3,i)
%    h=contour(X,Y,reshape(xv_ds(:,b30(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',[])
%    title(num2str(i))
% end
% suptitle('Original MSLP')
% 
% 
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 8 6]; % [left bottom width height]
% %print(['/home/nsteiger/Dropbox/easterisle/ei_underlyingdata_soms' num2str(somRows*somCols) '_mslp_0_' timestp2 '.png'],'-dpng','-r300');





% %--------------------------------------------------------------
% % What SOMs occur during big rain events vs small rain events?
% %--------------------------------------------------------------
% 
% % Frequency of BMUs during megadroughts
% sm_pfrq=histcounts(bmus(smrn),0.5:1:(somRows*somCols+0.5))/length(smrn);
% bg_pfrq=histcounts(bmus(bgrn),0.5:1:(somRows*somCols+0.5))/length(bgrn);
% 
% sm_dfq=(sm_pfrq-pfrq)./pfrq;
% disp('Small SOM p freq change')
% disp(sm_dfq)
% bg_dfq=(bg_pfrq-pfrq)./pfrq;
% disp('Big SOM p freq change')
% disp(bg_dfq)
% 
% 
% % Resample the BMUs to get uncertainty bounds for an "expected" range of frequencies
% rsz=500;
% s_rf=zeros(rsz,somRows*somCols);
% b_rf=zeros(rsz,somRows*somCols);
% for i=1:rsz
%    % Sample the reference period to get ranges for droughts and pluvials
%    srs=randsample(bmus,length(smrn),'true');
%    brs=randsample(bmus,length(bgrn),'true');
%    s_rf(i,:)=histcounts(srs,0.5:1:(somRows*somCols+0.5))/length(smrn);
%    b_rf(i,:)=histcounts(brs,0.5:1:(somRows*somCols+0.5))/length(bgrn);
% end
% 
% % Which values fall outside the expected 2 sigma range?
% x=1:somRows*somCols;
% s1=mean(s_rf)+2*std(s_rf);
% s2=mean(s_rf)-2*std(s_rf);
% sm_ux=find(sm_pfrq>s1 | sm_pfrq<s2);
% sm_in=setdiff(x,sm_ux);
% 
% b1=mean(b_rf)+2*std(b_rf);
% b2=mean(b_rf)-2*std(b_rf);
% bg_ux=find(bg_pfrq>b1 | bg_pfrq<b2);
% bg_in=setdiff(x,bg_ux);
% 
% 
% figure
% hold on
% % Frequency from distributions
% errorbar(1:somRows*somCols,mean(s_rf),2*std(s_rf),'o','markeredgecolor',[0.5 0.5 0.5],...
% 'color',[0.5 0.5 0.5],'linewidth',2)
% 
% % Frquency from actual with outliers differently colored
% plot(x(sm_in),sm_pfrq(sm_in),'p','markersize',10,'color',clrs.black,'markerfacecolor',clrs.black)
% plot(x(sm_ux),sm_pfrq(sm_ux),'p','markersize',10,'color',clrs.red,'markerfacecolor',clrs.red)
% 
% legend('Reference','Expected','Unexpected','location','eastoutside')
% xlabel('Node number','fontsize',18)
% ylabel('Frequency','fontsize',18)
% xticks(1:somRows*somCols)
% xtickangle(90)
% xlim([0.5 somRows*somCols+0.5])
% %ax = gca;ax.YAxisLocation='right';
% ylim([0 0.08])
% title('Small rainfall events')
% set(gca,'fontsize',16)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 8 5]; % [left bottom width height]
% box on
% 
% print('-dpng',['/home/nsteiger/Dropbox/easterisle/smallrain_som' num2str(somRows*somCols) '_ptrnfrq_' timestp2 '.png'],'-r300')
% 
% 
% 
% figure
% hold on
% % Frequency from distributions
% errorbar(1:somRows*somCols,mean(b_rf),2*std(b_rf),'o','markeredgecolor',[0.5 0.5 0.5],...
% 'color',[0.5 0.5 0.5],'linewidth',2)
% 
% % Frquency from actual with outliers differently colored
% plot(x(bg_in),bg_pfrq(bg_in),'p','markersize',10,'color',clrs.black,'markerfacecolor',clrs.black)
% plot(x(bg_ux),bg_pfrq(bg_ux),'p','markersize',10,'color',clrs.red,'markerfacecolor',clrs.red)
% 
% legend('Reference','Expected','Unexpected','location','eastoutside')
% xlabel('Node number','fontsize',18)
% ylabel('Frequency','fontsize',18)
% xticks(1:somRows*somCols)
% xtickangle(90)
% xlim([0.5 somRows*somCols+0.5])
% %ax = gca;ax.YAxisLocation='right';
% ylim([0 0.08])
% title('Large rainfall events')
% set(gca,'fontsize',16)
% fig = gcf;fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 2 8 5]; % [left bottom width height]
% box on
% 
% 
% print('-dpng',['/home/nsteiger/Dropbox/easterisle/largerain_som' num2str(somRows*somCols) '_ptrnfrq_' timestp2 '.png'],'-r300')
% 
% 
% close all
% save(['soms_analysis_estrisl_' timestp2 '.mat'],'-v7.3')














