%===========================================
% 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','xmeta')

clat=xmeta{1}.lat;
clon=xmeta{1}.lon;

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

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

% Which proxies to use?
%prxytyp='p'; % PDSI-like trees
prxytyp='a'; % all proxies

pd=[];j=1;
rv=[];
pcr=[];
plt=[];pln=[];
psmthd=cell(length(pi),1);
for i=1:length(pi)
   if strcmp(prxytyp,'p')
      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
   elseif strcmp(prxytyp,'a')
      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

prxyrs=proxy{1}.age_yng;

% LOAD OBSERVATIONAL DATA
%mon_avg_o=12;mon_avg_f=2;
%obs_dtp1='vanderS_pdsi';
mon_avg_o=4;mon_avg_f=3;
obs_dtp1='bearth'; % particular data type specified by input
[O_t,xlat,xlon,x_yrs] = load_obs_S(obs_dtp1,mon_avg_o,mon_avg_f);

% Get locations of proxies in obs
if ismember(obs_dtp1,{'bearth','hadcrut4'})
    [pl_o] = nearest_latlon(xlat,xlon,plt,pln);
else % mask out land or ocean
    mask=sum(~isnan(O_t),2);
    mask(mask==0)=NaN;
    mask(~isnan(mask))=1;
    [pl_o] = nearest_latlon(xlat,xlon,plt,pln,mask);
end

%analyrs1=1880:1980;
%[~,~,ir1] = intersect(analyrs1,r_o:r_f); % indices of recon segment

% Decadal avg
%O_t=ann2xt(O_t,10);

% Obs data at proxy locations
Opd=O_t(pl_o,:);
% Correlation matrix [lat*lon,proxy]
crm_f=corr(O_t',Opd','rows','pairwise');
% Max squared correlation
mx_crm=max(crm_f.^2,[],2);


load('./colormaps/MPL_BuPu.mat')
s=load('coast');mlat=s.lat;mlon=s.long;

%pi=find(p_lat<0 & p_lat>-57 & p_lon>270 & p_lon<320);
figure
%h=axesm('MapProjection','wagner4','MapLatLimit',[-60 10],'MapLonLimit',[270 330]);
%h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 15],'MapLonLimit',[270 330]);
h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 15],'MapLonLimit',[270 330],...
'grid','on','MLineLocation',[clon(1:72) (clon(73:end)-360)],'PLineLocation',clat,'GLineStyle','-');
%h=axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
pcolorm(xlat,xlon,reshape(mx_crm,length(xlat),length(xlon)))
%colormap(cmap)
caxis([0 1])
hold on
plotm(mlat,mlon,'color',[0 0 0])
%title(['Mean r^2 map at ' num2str(tsl(i)) ' CE'])
set(gca,'FontSize',14)
colorbar
tightmap
print(['/home/nsteiger/Dropbox/sa-phyda/phyda_proxy_maxcorr_map_t2m_allproxies_y' num2str(length(pl_o))],'-dpng','-r300')
%print(['/home/nsteiger/Dropbox/sa-phyda/phyda_proxy_maxcorr_map_pdsi_moistproxies_y' num2str(length(pl_o))],'-dpng','-r300')
%print(['/home/nsteiger/Dropbox/sa-phyda/phyda_proxy_maxcorr_decadal_map_pdsi_moistproxies_y' num2str(length(pl_o))],'-dpng','-r300')
























