
addpath('/home/nsteiger/MATLAB/matlab-goody-bag/cptcmap-pkg-master/cptcmap')
load distinct_colors_20.mat

%===========================================
% WHAT DO THE LOCAL SA PROXIES DO OVER TIME
%===========================================
efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:40:23.mat';
pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';
load([pth,efl],'proxy')

% Find proxies over southern South America
p_lat=zeros(length(proxy),1);
p_lon=zeros(length(proxy),1);
for i=1:length(proxy)
   p_lat(i)=proxy{i}.location(1);
   p_lon(i)=proxy{i}.location(2);
end
% Find which proxies are over south america land areas
%pi=find(p_lat<0 & p_lat>-57 & p_lon>270 & p_lon<320);
% Find which proxies are in the central Chile region
pi=find(p_lat<-26 & p_lat>-45 & p_lon>270 & p_lon<293);

% Pick other regions
%pi=find(p_lat<32 & p_lat>0 & p_lon>200 & p_lon<320); % mesoamerica
%pi=find(p_lat<90 & p_lat>32 & p_lon>200 & p_lon<320); % north america
%pi=find(p_lat<45 & p_lat>0 & p_lon>90 & p_lon<145); % monsoon asia
%pi=find(p_lat<72 & p_lat>28 & p_lon>0 & p_lon<45); % europe 1
%pi=find(p_lat<72 & p_lat>28 & p_lon>344 & p_lon<360); % europe 2


% Plot locations of proxies as a check
% s=load('coast');mlat=s.lat;mlon=s.long;
% 
% figure
% h=axesm('MapProjection','eqdcylin');
% hold on
% plotm(mlat,mlon,'color',[0.5 0.5 0.5])
% plotm(p_lat(pi),p_lon(pi),'^');

pd=[];j=1;
rv=[];
pcr=[];
plt=[];pln=[];
psmthd=cell(length(pi),1);
for i=1:length(pi)
   if strcmp(proxy{pi(i)}.ptype,'Tree Rings')
      if strcmp(proxy{pi(i)}.psm_mthd,'p')
      %if strcmp(proxy{pi(i)}.psm_mthd,'t')
	 pd(:,j)=proxy{pi(i)}.data;
	 rv(j)=proxy{pi(i)}.R;
	 pcr(j)=proxy{pi(i)}.pcr;
	 plt(j)=proxy{pi(i)}.location(1);
	 pln(j)=proxy{pi(i)}.location(2);
	 j=j+1;
      end
   end
end

prxyrs=proxy{1}.age_yng;


% %=======================
% % COMPUTE WEIGHTED MEAN
% %=======================
% 
% % Make sure the sign of the proxy data lines up with positive correlation with PDSI
% pds=bsxfun(@times,pd,sign(pcr));
% 
% %mean = sum(weights.*data)/sum(weights)
% %wt=1./rv;
% wt=abs(pcr); % weight according to proxy correlation with obs
% 
% enz=1000;
% pwmne=zeros(size(pd,1),enz);
% for i=1:enz
%    % Random sample of proxies
%    ri=randsample(size(pd,2),round(size(pd,2)*0.5));
%    % Compute the composite
%    rcn0=nansum(bsxfun(@times,wt(ri),pd(:,ri)),2)/sum(wt(ri));
%    % Scale the composite
% error('NOT FINISHED YET...')
%    %pwmne(:,i)=
% end
% 
% 
% figure
% plot(prctile(pwmne,[5],2),'color',clrs.grey,'linewidth',1.5)
% hold on
% plot(prctile(pwmne,[95],2),'color',clrs.grey,'linewidth',1.5)
% plot(prctile(pwmne,[50],2),'color',clrs.purple,'linewidth',2)
% xlim([1000 2000])
% set(gca,'fontsize',16)
% xlabel('Year','fontsize',18)
% ylabel('Weighted mean','fontsize',18)
% 
% %print(['/home/nsteiger/Dropbox/sa-phyda/sa_proxy_data_wmean.png'],'-dpng','-r300');

%==========================
% COMPUTE POST 1800 TRENDS
%==========================


% Make sure the sign of the proxy data lines up with positive correlation with PDSI
pds=bsxfun(@times,pd,sign(pcr));

% What years to look at?
analyrs1=1800:2000;
[~,~,ir1] = intersect(analyrs1,prxyrs); % indices of recon segment

% Make sure the sign of the proxy data lines up with positive correlation with PDSI
%pds=bsxfun(@times,pd,sign(pcr));

% Compute the trends
[trs]=trend(pds(ir1,:),prxyrs(ir1),'omitnan');



% DENSITY OF TRENDS

figure
histogram(trs,19)
hold on
yl=ylim;
plot(zeros(10,1),linspace(yl(1),yl(2),10),'k--','linewidth',2)
ylim(yl)
set(gca,'fontsize',16)
xlabel('Trend','fontsize',18)
ylabel('Count','fontsize',18)

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

fpth='/home/nsteiger/Dropbox/sa-phyda/';
%print([fpth 'sa_proxy_trend_pdf_18002000.png'],'-dpng','-r300');

%-------------------------
% DENSITY OF proxy values
%-------------------------

% look at 100 year means of each proxy time series
mpvals=zeros(length(4:9),size(pds,2));
k=1;
for i=4:9
   v1=[1 i 0 0];v1=sscanf(sprintf('%d',v1),'%li');
   v2=[1 i 9 9];v2=sscanf(sprintf('%d',v2),'%li');
   [~,~,ir] = intersect(v1:v2,prxyrs); % indices of recon segment
   for j=1:size(pds,2)
      mpvals(k,j)=nanmean(pds(ir,j));
   end
   k=k+1;
end

c0={'1400s','1500s','1600s','1700s','1800s','1900s'};

% Colors for boxplots
cm2=cat(1,clrs.grey,clrs.grey,clrs.grey,clrs.grey,clrs.grey,clrs.blue);


% distributions of the values
figure
hold on
%hl=boxplot(r2',c0','symbol','','colors',cm2);
%hl=boxplot(mpvals','Labels',c0,'symbol','','colors',cm2);
hl=boxplot(mpvals','Labels',c0,'symbol','o','colors',cm2);
set(hl,'linew',2)
set(findobj(gcf,'-regexp','Tag','\w*Whisker'),'LineStyle','-')% whisker lines solid
ax=gca;
ylim([-1 1])
%ax.YGrid = 'on';
line('xdata',ax.XLim,'ydata',[nanmedian(mpvals(end,:)) nanmedian(mpvals(end,:))],...
'linestyle','-','color',clrs.red,'linewidth',1)
line('xdata',ax.XLim,'ydata',[prctile(mpvals(end,:),25) prctile(mpvals(end,:),25)],...
'linestyle','--','color',clrs.black,'linewidth',1)
line('xdata',ax.XLim,'ydata',[prctile(mpvals(end,:),75) prctile(mpvals(end,:),75)],...
'linestyle','--','color',clrs.black,'linewidth',1)

ylabel('Proxy values (std)','fontsize',18)
set(gca,'FontSize',16,'XTickLabelRotation',45)
fig = gcf;fig.PaperUnits = 'inches';
%fig.PaperPosition = [0 2 4.5 4.5]; % [left bottom width height]
fig.PaperPosition = [0 2 3.5 3]; % add more vars

print(['/home/nsteiger/Dropbox/sa-phyda/sana_proxy_100yrmeans.png'],'-dpng','-r300');



return


% Density of values 
[fm,xim]=ksdensity(pvls_mdrn(:));
[f6,xi6]=ksdensity(pvls_1600(:));
[f4,xi4]=ksdensity(pvls_1400(:));

figure
hold on
plot(xim,fm,'linewidth',2,'color',clrs.black)
plot(xi6,f6,'linewidth',2,'color',clrs.orange)
plot(xi4,f4,'linewidth',2,'color',clrs.blue)
set(gca,'fontsize',16)
xlabel('Prxy vals (std units)','fontsize',18)
ylabel('Density','fontsize',18)
%legend('1900-2000','1600-1800','1400-1599','location','eastoutside')
box on
grid on
xlim([-5 5])

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

print([fpth 'sa_proxyvals_comp_pdf.png'],'-dpng','-r300');


return



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

figure
hold on
%C=importdata('redblue_256_rgb.txt');
%colormap(flip(C,1)/255)
%cptcmap('GMT_drywet','flip',true)
cptcmap('GMT_red2green')
caxis([-0.007 0.007])
%set(get(colorbar('location','southoutside'),'xlabel'),'string',...
%    'Trend per year','fontsize',14)
colorbar
set(gca,'fontsize',14)
axesm('MapProjection','wagner4','MapLatLimit',[-60,-15],'MapLonLimit',[-83,-37]);
%title('VSL Pseudotrees','fontsize',18)
plotm(mlat,mlon,'color',[0.5 0.5 0.5])
hs1=scatterm(plt,pln,30,trs,'filled');
%hs2=scatterm(tlat(sh),tlon(sh),15,mL_djf(sh),'filled');
set(get(hs1,'Children'),'MarkerEdgeColor','k');
%set(get(hs2,'Children'),'MarkerEdgeColor','k');
%alpha(0.1);
tightmap
hold off
title('Trend 1800-2000 PDSI-like trees')

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

fpth='/home/nsteiger/Dropbox/sa-phyda/';
%print([fpth 'sa_proxy_trend_map_18002000.png'],'-dpng','-r300');


% DENSITY OF TRENDS
[f,xi]=ksdensity(trs);

figure
plot(xi,f,'linewidth',3)
hold on
yl=ylim;
plot(zeros(10,1),linspace(yl(1),yl(2),10),'k--','linewidth',2)
ylim(yl)
set(gca,'fontsize',16)
xlabel('Trend (std units per year)','fontsize',18)
ylabel('Density','fontsize',18)

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

%print([fpth 'sa_proxy_trend_pdf_18002000.png'],'-dpng','-r300');









