
% Compute the spatial skill

tic

%addpath('../')

% ALL PROXY TYPES, LMR DATABASE
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_AprMar_prxydtst_5_prxtp_tca_2604_swtchbld100_21-Dec-2017_16:05:26.mat';
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_tca_2598_swtchbld100_21-Dec-2017_16:19:04.mat';
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2588_swtchbld100_21-Dec-2017_16:30:40.mat';

%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:33:07.mat';
efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:40:23.mat';
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:49:45.mat';

disp(['Loading data: ' efl])

pth='./output-recon/'; load([pth,efl])

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

%load(['/d1/nsteiger/climate-data/cesm-lme/cesm_lme_010_pdsi_output_biascorr_AWC_c_7_u2_gcm.mat'],'lndidx')

k=find(state_tp=='q'); % PDSI
Xa_pdsi_lndonly=Xa_m(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);
lndidx=xmeta{k}.lndidx;
lat=xmeta{k}.lat;lon=xmeta{k}.lon;
Xa_pdsi_f=NaN(length(lat)*length(lon),reconYrs);
Xa_pdsi_f(lndidx,:)=Xa_pdsi_lndonly;
% Additional fields for CRPSS
Xa_pdsi_s=Xa_sigma(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:); % pull out sigma for CRPSS
Xa_pdsi_sigf=NaN(length(lat)*length(lon),reconYrs);
Xa_pdsi_sigf(lndidx,:)=Xa_pdsi_s;
Xb_pdsi_0=Xb_o(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:); % pull out sigma for CRPSS
Xb_pdsi=NaN(length(lat)*length(lon),length(p_yrs));
Xb_pdsi(lndidx,:)=Xb_pdsi_0;

k=find(state_tp=='j'); % SPEI
Xa_spei_lndonly=Xa_m(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);
Xa_spei_f=NaN(length(lat)*length(lon),reconYrs);
Xa_spei_f(lndidx,:)=Xa_spei_lndonly;
% Additional fields for CRPSS
Xa_spei_s=Xa_sigma(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:); % pull out sigma for CRPSS
Xa_spei_sigf=NaN(length(lat)*length(lon),reconYrs);
Xa_spei_sigf(lndidx,:)=Xa_spei_s;
Xb_spei_0=Xb_o(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:); % pull out sigma for CRPSS
Xb_spei=NaN(length(lat)*length(lon),length(p_yrs));
Xb_spei(lndidx,:)=Xb_spei_0;



k=find(state_tp=='t'); % T2m
Xa_t2m=Xa_m(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:);
% Additional fields for CRPSS
Xa_t2m_sig=Xa_sigma(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:); % pull out sigma for CRPSS
Xb_t2m=Xb_o(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:); % pull out sigma for CRPSS





% Analysis options
anal_da='n'; % drought atlas verification?
anal_obs='y'; % observational verification?
plot_p='n'; % plot proxies?
make_figs='n'; % Save figures?


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

%=========================================
% VERIFICATION AGAINST OBSERVATIONAL DATA
%=========================================

sptlskll=cell(0);



%====================
% T 2m  VERIFICATION
%====================

disp('Verifying against T2m...')
toc

% Load PDSI
[Xobs,xlat,xlon,x_yrs] = load_obs_S('bearth',mon_avg_o,mon_avg_f);
da=reshape(Xobs,length(xlat),length(xlon),length(x_yrs));

% Interpolate to CESM resolution
%ndm=load('dai_pdsi_mask.mat');
% create grids for interpolation
[xn,yn]=meshgrid(xlon,xlat);
[xc,yc]=meshgrid(lon,lat);
% note masks of both data types
%nm=find(~isnan(da(:,:,end)));
%cm=find(ndm.sftlf==1);
%da(da<-99)=NaN; % put nans in fill spots but after determining mask

disp('Interplating obs to climate model res...')
da_i=nan(length(lat),length(lon),size(da,3));
for i=1:size(da,3)
   % linear interpolation and extrapolation for irregular data
   da1=da(:,:,i);
   F = scatteredInterpolant(xn(:),yn(:),da1(:),'linear','linear');
   vq=F(xc,yc); % do the interpolation
   da_i(:,:,i)=vq;
end

% COMPUTE SKILL OF THE RECONSTRUCTION

% get years of overlap
%[~,ir,ih] = intersect(r_o:r_f,time);
% Only look at a certain segment of the reconstructions
%analyrs=1800:1899;
analyrs=1870:1990;
[~,~,ir] = intersect(analyrs,r_o:r_f); % indices of recon segment
[~,~,id] = intersect(analyrs,x_yrs); % indices of DA segment

% Anomaly with respect to analysis period
dai=reshape(da_i,length(lat)*length(lon),size(da,3));

ph_nm=zeros(size(Xa_t2m,1),length(ir)); 
be_nm=zeros(size(dai,1),length(ir));
for j=1:length(lat)*length(lon)
   % Remove mean
   ph_nm(j,:)=Xa_t2m(j,ir)-nanmean(Xa_t2m(j,ir));
   be_nm(j,:)=dai(j,id)-nanmean(dai(j,id));
end

% Just look at Eurasia
pn=find(lon>0 & lon<150);
pt=find(lat>40 & lat<70);

p0=reshape(ph_nm,length(lat),length(lon),length(ir));
b0=reshape(be_nm,length(lat),length(lon),length(ir));
ph_nm1=reshape(p0(pt,pn,:),length(pn)*length(pt),length(ir));
be_nm1=reshape(b0(pt,pn,:),length(pn)*length(pt),length(ir));

% PLOT MEAN TIME SERIES AGAINST EACH OTHER
figure
hold on
plot(ir,wmean(p0(pt,pn,:),lat(pt)),'-or','linewidth',2)
plot(ir,wmean(b0(pt,pn,:),lat(pt)),'-ok','linewidth',2)
ax=gca;
line('ydata',ax.YLim,'xdata',[1883 1883],...
'linestyle','-','color','green','linewidth',2)
ax.YGrid = 'on';
box on
xlim([1871 1891])
ylim([-3 3])
xlabel('Year','fontsize',16)
ylabel('Anomaly correlation','fontsize',16)
set(gca,'FontSize',14)
legend('PHYDA','BEarth')
title('DJF Eurasia')

%print(['/home/nsteiger/Dropbox/enso-volcanoes/phyda_bearth_djf_mean_18711891_eurasia.png'],'-dpng','-r300');



% COMPUTE ANOMALY CORRELATION TO SEE IF DATING IS CORRECT
spt_corr0=zeros(length(ir)-2,1);
spt_corrm1=zeros(length(ir)-2,1);
spt_corrp1=zeros(length(ir)-2,1);
k=1;
for i=2:(length(ir)-1)
   spt_corr0(k)=nancorr(ph_nm1(:,i),be_nm1(:,i));
   spt_corrm1(k)=nancorr(ph_nm1(:,i-1),be_nm1(:,i));
   spt_corrp1(k)=nancorr(ph_nm1(:,i+1),be_nm1(:,i));
   k=k+1;
end

figure
hold on
plot(ir(2:end-1),spt_corr0,'-ok','linewidth',2)
plot(ir(2:end-1),spt_corrm1,'-ob','linewidth',2)
plot(ir(2:end-1),spt_corrp1,'-or','linewidth',2)
ax=gca;
line('ydata',ax.YLim,'xdata',[1883 1883],...
'linestyle','-','color','green','linewidth',2)
ax.YGrid = 'on';
box on
xlim([1871 1891])
xlabel('Year','fontsize',16)
ylabel('Anomaly correlation','fontsize',16)
set(gca,'FontSize',14)
legend('0','-1','+1')
title('AnomCorr(PHYDA,BEarth) Eurasia')

%print(['/home/nsteiger/Dropbox/enso-volcanoes/phyda_bearth_djf_anomaly_corr_18711891_eurasia.png'],'-dpng','-r300');


% COMPUTE CRPS TO SEE IF DATING IS CORRECT
p0=reshape(Xa_t2m_sig(:,ir),length(lat),length(lon),length(ir));
phsg=reshape(p0(pt,pn,:),length(pn)*length(pt),length(ir));

spt_crps0=zeros(length(ir)-2,1);
spt_crpsm1=zeros(length(ir)-2,1);
spt_crpsp1=zeros(length(ir)-2,1);
k=1;
for i=2:(length(ir)-1)
   spt_crps0(k)=crps_m(ph_nm1(:,i),phsg(:,i),be_nm1(:,i));
   spt_crpsm1(k)=crps_m(ph_nm1(:,i-1),phsg(:,i-1),be_nm1(:,i));
   spt_crpsp1(k)=crps_m(ph_nm1(:,i+1),phsg(:,i+1),be_nm1(:,i));
   k=k+1;
end

figure
hold on
plot(ir(2:end-1),spt_crps0,'-ok','linewidth',2)
plot(ir(2:end-1),spt_crpsm1,'-ob','linewidth',2)
plot(ir(2:end-1),spt_crpsp1,'-or','linewidth',2)
ax=gca;
line('ydata',ax.YLim,'xdata',[1883 1883],...
'linestyle','-','color','green','linewidth',2)
ax.YGrid = 'on';
box on
xlim([1871 1891])
ylim([0.5 3])
xlabel('Year','fontsize',16)
ylabel('Spatial CRPS','fontsize',16)
set(gca,'FontSize',14)
legend('0','-1','+1')
title('CRPS(PHYDA,BEarth) Eurasia')

print(['/home/nsteiger/Dropbox/enso-volcanoes/phyda_bearth_djf_crps_18711891_eurasia.png'],'-dpng','-r300');









% % Just look at central Pacific
% pt=find(lat<25 & lat>-25);
% pn=find(lon>130 & lon<280); % corresponds to 130 E to 70 W
% 
% 
% p0=reshape(ph_nm,length(lat),length(lon),length(ir));
% b0=reshape(be_nm,length(lat),length(lon),length(ir));
% ph_nm1=reshape(p0(pt,pn,:),length(pn)*length(pt),length(ir));
% be_nm1=reshape(b0(pt,pn,:),length(pn)*length(pt),length(ir));
% 
% 
% % COMPUTE ANOMALY CORRELATION TO SEE IF DATING IS CORRECT
% spt_corr0=zeros(length(ir)-2,1);
% spt_corrm1=zeros(length(ir)-2,1);
% spt_corrp1=zeros(length(ir)-2,1);
% k=1;
% for i=2:(length(ir)-1)
%    spt_corr0(k)=nancorr(ph_nm1(:,i),be_nm1(:,i));
%    spt_corrm1(k)=nancorr(ph_nm1(:,i-1),be_nm1(:,i));
%    spt_corrp1(k)=nancorr(ph_nm1(:,i+1),be_nm1(:,i));
%    k=k+1;
% end
% 
% figure
% hold on
% plot(ir(2:end-1),spt_corr0,'-ok','linewidth',2)
% plot(ir(2:end-1),spt_corrm1,'-ob','linewidth',2)
% plot(ir(2:end-1),spt_corrp1,'-or','linewidth',2)
% ax=gca;
% line('ydata',ax.YLim,'xdata',[1883 1883],...
% 'linestyle','-','color','green','linewidth',2)
% ax.YGrid = 'on';
% box on
% xlim([1871 1891])
% xlabel('Year','fontsize',16)
% ylabel('Anomaly correlation','fontsize',16)
% set(gca,'FontSize',14)
% legend('0','-1','+1')
% title('AnomCorr(PHYDA,BEarth) central Pacific')
% 
% print(['/home/nsteiger/Dropbox/enso-volcanoes/phyda_bearth_djf_anomaly_corr_18711891_pacific.png'],'-dpng','-r300');
% 

return

% PLOT PARTICULAR YEARS


ph_nm=reshape(ph_nm,length(lat),length(lon),length(ir));
be_nm=reshape(be_nm,length(lat),length(lon),length(ir));


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

for j=1:length(ir)

   figure
   hold on
   %h=axesm('MapProjection','wagner4','MapLatLimit',[(min(xlat)) (max(xlat))],'MapLonLimit',[min(xlon) max(xlon)]);
   h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 max(lat)]);
   pcolorm(lat,lon,ph_nm(:,:,j)); % add additional value to avoid line in pcolor plots
   plotm(mlat,mlon,'color',[0.5 0.5 0.5])
   title(['PHYDA ' num2str(ir(j))])
   set(gca,'fontsize',18)
   %set(get(colorbar('location','southoutside'),'xlabel'),'string',sn{j},'fontsize',14)
   caxis([-3 3])
   %set(cbh,'YTick',[-1:0.2:1])
   %delete( cbh )
   hold off

   figure
   hold on
   %h=axesm('MapProjection','wagner4','MapLatLimit',[(min(xlat)) (max(xlat))],'MapLonLimit',[min(xlon) max(xlon)]);
   h=axesm('MapProjection','eqdcylin','MapLatLimit',[-60 max(lat)]);
   pcolorm(lat,lon,be_nm(:,:,j)); % add additional value to avoid line in pcolor plots
   plotm(mlat,mlon,'color',[0.5 0.5 0.5])
   title(['BEarth ' num2str(ir(j))])
   set(gca,'fontsize',18)
   %set(get(colorbar('location','southoutside'),'xlabel'),'string',sn{j},'fontsize',14)
   caxis([-3 3])
   %set(cbh,'YTick',[-1:0.2:1])
   %delete( cbh )
   hold off
end





return


disp('Computing CRPSS...')

% compute mean crpss
frcst_c=zeros(length(lat)*length(lon),1);
m_crpss=zeros(length(lat)*length(lon),1);
for j=1:length(lat)*length(lon)
   % Remove mean
   f_mu=Xa_t2m(j,ir)-nanmean(Xa_t2m(j,ir));
   y_o=dai(j,id)-nanmean(dai(j,id));
   frcst_c(j)=crps_m(f_mu,Xa_t2m_sig(j,ir),y_o);
 
   % Randomly draw from prior to get a mean CRPS
   ref_c0=zeros(250,1);
   parfor i=1:length(ref_c0)    
      rsxb=randsample(Xb_t2m(j,:),length(y_o));
      rsxb=rsxb-mean(rsxb);
      ref_c0(i)=crps_m(rsxb,std(rsxb),y_o);
   end
   ref_c=nanmean(ref_c0);

   % compute skill score
   m_crpss(j)=1-frcst_c(j)/ref_c;

end
% Save skill maps
sptlskll.t2m.rp=rp;
sptlskll.t2m.ce=ce;
sptlskll.t2m.crps=frcst_c;
sptlskll.t2m.crpss=m_crpss;



% Save the skill 
save(['./output-recon/spatial_skill_' efl],'-v7.3')





