
clear

addpath('../')

%==============================
% LOAD THE RECONSTRUCTION DATA
%==============================

% ALL PROXY TYPES, LMR DATABASE
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:49:45.mat';
efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:33:07.mat';
% INDEX ONLY RECONSTRUCTION WITH FULL ENSEMBLES
%efl='cesm_lme010_r12000_p998_state_2ozn_avg_AprMar_prxydtst_5_prxtp_tca_2605_swtchbld100_16-Nov-2017_15:57:59.mat';

%pth='/d2/nsteiger/output-da/hydroclimate/'; load([pth,efl])
pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';load([pth,efl])

% Temperature
k=find(state_tp=='t');
Xv=Xb_o(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);

lat=xmeta{k}.lat;
lon=xmeta{k}.lon;

% Time slices to analyze
tsl=[1000,1400,1800];

%======================================================
% COMPUTE THE CORRELATION MAPS FOR EACH PROXY LOCATION
%======================================================
plo=zeros(length(proxy),1);
pdata=NaN(length(proxy),3);
for i=1:length(proxy)
   plat=proxy{i}.location(1);
   plon=proxy{i}.location(2);
   plo(i)=nearest_latlon(lat,lon,plat,plon);
   % Proxy data at time slices
   pdata(i,:)=proxy{i}.data(tsl);
end
% Model data at proxy locations
Xvp=Xv(plo,:);
% Correlation matrix [lat*lon,proxy]
crm_f=corr(Xv',Xvp');

%-------------------------------------------------------------------
% WHAT DO THE CORRELATION MAPS LOOK LIKE AT DIFFERENT TIME PERIODS?
%-------------------------------------------------------------------

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

for i=1:length(tsl)
   
   prm=find(~isnan(pdata(:,i)));

   % Model data at proxy locations
   Xvp=Xv(plo(prm),:);
   % Correlation matrix [lat*lon,proxy]
   crm=corr(Xv',Xvp');
   % Mean squared correlation
   m_crm=mean(crm.^2,2);

   figure
   %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)]);
   pcolorm(lat,lon,reshape(m_crm,length(lat),length(lon)))
   colormap(cmap)
   caxis([0 0.15])
   hold on
   plotm(mlat,mlon,'color',[0.25 0.25 0.25])
   title(['Mean r^2 map at ' num2str(tsl(i)) ' CE'])
   colorbar
   tightmap
   print(['phyda_proxy_meancorr_map_y' num2str(tsl(i))],'-dpng','-r300')


end


for i=1:length(tsl)
   
   prm=find(~isnan(pdata(:,i)));

   % Model data at proxy locations
   Xvp=Xv(plo(prm),:);
   % Correlation matrix [lat*lon,proxy]
   crm=corr(Xv',Xvp');
   % Mean squared correlation
   m_crm=max(crm.^2,[],2);

   figure
   %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)]);
   pcolorm(lat,lon,reshape(m_crm,length(lat),length(lon)))
   colormap(cmap)
   caxis([0 1])
   hold on
   plotm(mlat,mlon,'color',[0.25 0.25 0.25])
   title(['Max r^2 map at ' num2str(tsl(i)) ' CE'])
   colorbar
   tightmap
   print(['phyda_proxy_maxcorr_map_y' num2str(tsl(i))],'-dpng','-r300')


end









