
% Different options when computing SOMs

% LAND + OCEAN
%fln='output_som4_wt_y_d_1_18-Apr-2018_16:17:50.mat';
%fln='output_som4_wt_y_d_2_18-Apr-2018_15:58:50.mat';
%fln='output_som6_wt_y_d_1_18-Apr-2018_15:29:20.mat';
%fln='output_som6_wt_y_d_2_18-Apr-2018_15:56:38.mat';
%fln='output_som8_wt_y_d_1_18-Apr-2018_17:03:04.mat';
%fln='output_som8_wt_y_d_2_18-Apr-2018_17:11:25.mat';
%fln='output_som9_wt_y_d_1_18-Apr-2018_16:54:38.mat';
%fln='output_som9_wt_y_d_2_18-Apr-2018_16:51:02.mat';
%fln='output_som12_wt_y_d_1_18-Apr-2018_17:22:36.mat';
%fln='output_som12_wt_y_d_2_18-Apr-2018_17:18:55.mat';
%fln='output_som12_wt_y_d_1_nodetrend_19-Apr-2018_14:25:02.mat';

% OCEAN ONLY
%fln='output_som4_d_1_01-May-2018_13:47:27.mat';
%fln='output_som6_d_1_02-May-2018_20:06:26.mat';
%fln='output_som8_d_1_02-May-2018_15:03:05.mat';
%fln='output_som9_d_1_01-May-2018_13:39:50.mat';
%fln='output_som12_d_1_02-May-2018_15:05:12.mat';

% For Hun
%fln='output_som12_25-Jul-2018_05:12:34.mat';


%fln0=fln;
%load(fln)

%---------------------------
%      PLOT PATTERNS
%---------------------------

load('./colormaps/distinct_colors_20.mat')
%load('./colormaps/WhiteYellowOrangeRed.mat')
C=load('./colormaps/NCV_blue_red.mat');cmap_br=C.cmap;
C=load('./colormaps/MPL_BrBG.mat');cmap_bg=C.cmap;


make_figs='n';
plot_sp='y';% make spatial plots?


analyrs4=800:1925;
[~,~,ir4] = intersect(analyrs4,r_o:r_f);

% PULL OUT SPATIAL MEAN OVER AMERICAN SOUTHWEST
%load('/d1/nsteiger/climate-data/cesm-lme/Nrth_A_Swthwst_mask.mat')
varpth=['/d6/haibo/lme/cam5/BLMTRC5CN/ALLFORCEINGS/010'];
sftlf=permute(ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'landmask'),[2 1]);
%lat=ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'lat');
%lon=ncread([varpth '/QSOIL/b.e11.BLMTRC5CN.f19_g16.010.clm2.h0.QSOIL.085001-184912.nc'],'lon');
% NASW = box of 235:255 E (125:105 W); 31:42 N (just the US, not Mexico)
% Find nearest lat/lon equal to or outside bounds
lonext=[235 255];latext=[31 42];
ln_a=find((abs(lon-lonext(1))-min(abs(lon-lonext(1))))<1e-5);
ln_b=find((abs(lon-lonext(2))-min(abs(lon-lonext(2))))<1e-5);
lt_a=find((abs(lat-latext(1))-min(abs(lat-latext(1))))<1e-5);
lt_b=find((abs(lat-latext(2))-min(abs(lat-latext(2))))<1e-5);
nasw=zeros(size(sftlf));
nasw(lt_a:lt_b,ln_a:ln_b)=sftlf(lt_a:lt_b,ln_a:ln_b);





if plot_sp=='y'


for i=1:somRows*somCols


   s=load('coast');mlat=s.lat;mlon=s.long;

   figure
   hold on
   %h=axesm('MapProjection','wagner4','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
   h1=axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
   %h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 max(lat)]);
   title(['Pattern ' num2str(i)],'fontsize',18,'Interpreter','none')
   pcolorm(lat,lon,reshape(ptrns(i,:),length(lat),length(lon)));
   %pcolorm(lat,[lon; 360],varZ); % add additional value to avoid line in pcolor plots
   plotm(mlat,mlon,'color',[0.5 0.5 0.5])
   colormap(h1,cmap_br) 
   set(gca,'fontsize',18)
   %set(get(colorbar('location','southoutside'),'xlabel'),'string',sn{j},'fontsize',14)
   %set(get(colorbar('location','southoutside'),'xlabel'),'string','std units','fontsize',18)
   caxis(h1,[-1.5 1.5])
   %cbh = findobj( 0, 'tag', 'Colorbar' ); % and
   %set(cbh,'YTick',[-1.5:0.25:1.5])
   %delete( cbh )
   tightmap
   hold off


end

end


 
 













% ORIGINAL VERSION OF THE CODE
% 
% if plot_p=='y'
%    % Pull out proxy locations
%    %plat=nan(length(proxy),1);
%    %plon=nan(length(proxy),1);
%    j=1;k=1;l=1;
%    for i=1:length(proxy)
%       if ismember(proxy{i}.ptype,{'tree','Tree Rings'})
% 	 tlat(j)=proxy{i}.location(1);
% 	 tlon(j)=proxy{i}.location(2);
% 	 j=j+1;
%       elseif ismember(proxy{i}.ptype,{'coral','sclerosponge','Corals and Sclerosponges','Bivalve'})
% 	 clat(k)=proxy{i}.location(1);
% 	 clon(k)=proxy{i}.location(2);
% 	 k=k+1;
%       else % Other proxy type 
% 	 olat(l)=proxy{i}.location(1);
% 	 olon(l)=proxy{i}.location(2);
% 	 l=l+1;
%       end 
%    end
% end
% 
% if plot_sp=='y'
% 
% for i=1:somRows*somCols
% 
% 
%    s=load('coast');mlat=s.lat;mlon=s.long;
% 
%    figure
%    hold on
%    %h=axesm('MapProjection','wagner4','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
%    h=axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
%    %h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 max(lat)]);
%    title(['Pattern ' num2str(i)],'fontsize',18,'Interpreter','none')
%    pcolorm(lat,lon,reshape(ptrns(i,:),length(lat),length(lon)));
%    %pcolorm(lat,[lon; 360],varZ); % add additional value to avoid line in pcolor plots
%    plotm(mlat,mlon,'color',[0.5 0.5 0.5])
%    if plot_p=='y';
%       h1=plotm(tlat,tlon,'s','color',clrs.green);
%       h2=plotm(clat,clon,'*','color',clrs.purple);
%       h3=plotm(olat,olon,'^','color',clrs.black);
%       %legend([h1,h2,h3],{'Tree-Rings','Corals','Others'},'location','southeast')
%    end 
%    colormap(cmap) % colormap from 'WhiteYellowOrangeRed'
%    set(gca,'fontsize',18)
%    %set(get(colorbar('location','southoutside'),'xlabel'),'string',sn{j},'fontsize',14)
%    set(get(colorbar('location','southoutside'),'xlabel'),'string','std units','fontsize',18)
%    caxis([-1.5 1.5])
%    %cbh = findobj( 0, 'tag', 'Colorbar' ); % and
%    %set(cbh,'YTick',[-1.5:0.25:1.5])
%    %delete( cbh )
%    hold off
% 
%    if  make_figs=='y'
%       print('-dpng',['./figs/som' num2str(somRows*somCols) '_ptrn_' num2str(i) '_' timestp2 '.png'],'-r300')
%    end
% 
% end
% 
% end
% 
% 
% % figure
% % for i=1:somRows*somCols
% %   subplot(somRows,somCols,i)
% %   h=pcolor(reshape(ptrns(i,:),length(lat),length(lon)));
% %   set(h,'edgecolor','none')
% % end
% % %ht=suptitle(['' fn(1:enddd-13)]);
% % %set(ht,'interpreter','none')
%  
%  
% %----------------------------------------------
% %  COMPARE SOM FREQUENCIES BY PERIOD
% %----------------------------------------------
% 
% 
% % Get BMUs for Mdroughts during 800-1600, for non-drought yrs 800-1600, all years 1601-1925
% 
% % BMU frequency during megadroughts (note droughts are determined analysis period relative)
% df=[];
% for i=1:length(ds0)
%   df0=ds0(i,1):ds0(i,2);
%   df=cat(2,df,df0);
% end
% dfc=histcounts(bmus(df),0.5:1:(somRows*somCols+0.5))/length(df);
% 
% % Frequency during non-droughts
% [~,ay1,~] = intersect(analyrs0,analyrs1);
% wf=setdiff(ay1(:),df(:));
% wfc=histcounts(bmus(wf),0.5:1:(somRows*somCols+0.5))/length(wf);
% 
% % Frequency during 1600-1925 period
% [~,ay2,~] = intersect(analyrs0,analyrs2);
% a2fc=histcounts(bmus(ay2),0.5:1:(somRows*somCols+0.5))/length(ay2);
%  
% % Frequency for each drought
% ds_pfrq=zeros(length(ds0),somRows*somCols);
% for i=1:length(ds0)
%   ds_pfrq(i,:)=histcounts(bmus(ds0(i,1):ds0(i,2)),0.5:1:(somRows*somCols+0.5))/length(ds0(i,1):ds0(i,2));
% end
% 
% disp(fln0)
% disp('Drought percent change of the frequency of the pattern')
% pd=(dfc-a2fc)./a2fc;
% disp(pd)
% 
% %----------------------------
% % PLOT BMUS FOR EACH PATTERN
% %----------------------------
% 
% C=load('./colormaps/MPL_BrBG.mat');
% 
% figure
% hold on
% % Mean freq during mega D
% b=bar(1:somRows*somCols,a2fc);
% b.FaceColor=[0.5 0.5 0.5];
% % Mean freq 1601-1925
% plot(1:somRows*somCols,dfc,'o','linewidth',2,'color',C.cmap(1,:),'markerfacecolor',C.cmap(1,:))
% plot(1:somRows*somCols,wfc,'^','linewidth',2,'color',C.cmap(end,:))
% legend('1601-1925','Megadrought 800-1600','Non-drought 800-1600','location','south')
% %title('Ptrn. Freq. MDrought')
% xlabel('Pattern number','fontsize',18)
% ylabel('Frequency','fontsize',18)
% xticks(1:somRows*somCols)
% set(gca,'fontsize',16)
% box on
% if  make_figs=='y'
%    print('-dpng',['./figs/som' num2str(somRows*somCols) '_ptrnfrq_' timestp2 '.png'],'-r300')
% end
% 
% 
% 
% %figur
% %plot(analyrs0,bmus,'or')
% %ylabel('Pattern')
% %yticks(1:6)
% %xlabel('Year')
% %box on
% %if  make_figs=='y'
% %   print('-dpng',['./figs/ptrn_frq_timeseries' timestp2 '.png',])
% %end
% 
% % Model the transition frequencies as Markov chain
% 
% % Compute transition matrix
% %[T1] = transM(bmus,1);
% 
% % Markov chain object
% %mc = dtmc(T1);
% 
% % Directed graph with color = probabilities
% %figure
% %graphplot(mc,'ColorEdges',true);
% 
% 
% 
% %----------------------------------------------
% %  COMPOSITES OF PDSI OVER THE SOMS
% %----------------------------------------------
% 
% if plot_sp=='y'
% 
% % Detrend data for composites 
% D0=reshape(Xpdsi(:,:,ir0),length(lat)*length(lon),length(ir0));
% dt_pdsi=detrend(D0')';% Detrend the data, if desired
% 
% 
% for i=1:somRows*somCols
% 
%    cind=find(bmus==i);
%    varX=reshape(mean(dt_pdsi(:,cind),2),length(lat),length(lon));
% 
%    figure
%    hold on
%    s=load('coast');mlat=s.lat;mlon=s.long;
%    %h=axesm('MapProjection','wagner4','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
%    %h=axesm('MapProjection','wagner4','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[120 33]);
%    h=axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
%    varZ=cat(2,varX,varX(:,1));
%    pcolorm(lat,[lon; 360],varZ);
%    title(['Pattern ' num2str(i)],'fontsize',18,'Interpreter','none')
%    %cptcmap('GMT_red2green')
%    %cptcmap('GMT_drywet')
%    %C=importdata('MPL_BrBG.txt');
%    C=load('./colormaps/MPL_BrBG.mat');
%    colormap(C.cmap)
%    caxis([-1.25 1.25])
%    set(gca,'fontsize',16)
%    set(get(colorbar('location','southoutside'),'xlabel'),'string','PDSI','fontsize',16)
%    cbh = findobj( 0, 'tag', 'Colorbar' );
%    %set(cbh,'YTick',[-1:0.2:1])
%    plotm(mlat,mlon,'color',[0.5 0.5 0.5])
%    tightmap
% 
%    if make_figs=='y';print('-dpng',['./figs/cmpst_pdsi_som' num2str(somRows*somCols) '_ptrn_' num2str(i) '_' timestp2 '.png'],'-r300'); end
% end
% 
% end
% 







