
% Compute the spatial skill


addpath('../')

% ALL PROXY TYPES, LMR DATABASE
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_AprMar_prxydtst_5_prxtp_tca_2605_swtchbld100_16-Nov-2017_16:03:43.mat';
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_DecFeb_prxydtst_5_prxtp_tca_2601_swtchbld100_16-Nov-2017_16:11:22.mat';
%efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_JunAug_prxydtst_5_prxtp_tca_2590_swtchbld100_16-Nov-2017_16:20:57.mat';

efl='cesm_lme010_r12000_p998_state_t2qjozn_avg_AprMar_prxydtst_5_prxtp_tca_2896_swtchbld100_05-Jan-2018_14:33:07.mat';

disp(['Loading data: ' efl])

pth='/d2/nsteiger/output-da/hydroclimate/subrecons_5Jan18/';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=='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
lat=xmeta{k}.lat;
lon=xmeta{k}.lon;

k=find(state_tp=='q');
Xvar=squeeze(Xa_m(xmeta{k}.id_X(1):xmeta{k}.id_X(2),:));
Xpdsi=NaN(length(xmeta{k}.lat)*length(xmeta{k}.lon),reconYrs);
Xpdsi(xmeta{k}.lndidx,:)=Xvar;
%Xpdsi=reshape(Xpdsi,length(xmeta{k}.lat),length(xmeta{k}.lon),reconYrs);




lndidx=xmeta{3}.lndidx;
ocnidx = setdiff((1:length(lat)*length(lon))',lndidx);
W=cosd(repmat(lat,[1 length(lon)]));%W=W(:);

% Analysis options
prnt_fgs='n'; % Save figures?

analyrs1=1000:2000;
analyrs2=1950:2000;
[~,~,ir1] = intersect(analyrs1,r_o:r_f); % indices of recon segment
[~,~,ir2] = intersect(analyrs2,r_o:r_f); % indices of recon segment


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

[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('Interpolating 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

dai=reshape(da_i,length(lat)*length(lon),size(da,3));

% Pull out just SSTs
d_sst=dai;
d_sst(lndidx,:)=NaN;
x_sst=Xa_t2m;
x_sst(lndidx,:)=NaN;

% Indices of the obs
[~,~,ib1] = intersect(analyrs1,x_yrs); 
[~,~,ib2] = intersect(analyrs2,x_yrs); 



%-------------
%  Load PDSI
%-------------

[Xobs,xlat,xlon,x_yrs] = load_obs_S('dai_pdsi',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(nm),yn(nm),da1(nm),'linear','linear');
   vq=F(xc(cm),yc(cm)); % do the interpolation
   % expand out the interpolated data points onto the climate model mask
   vd=nan(length(lat),length(lon));
   vd(cm)=vq; % full interpolated field
   da_i(:,:,i)=vd;
end

d_pdsi=reshape(da_i,length(lat)*length(lon),size(da,3));

% Indices of the obs
[~,~,id1] = intersect(analyrs1,x_yrs); 
[~,~,id2] = intersect(analyrs2,x_yrs); 


% Calculate sliding decadal differences, see which are most highly correlated,
% and see which decades give the largest decadal swings in temperatures 

% compute decadal moving average
%ws = 10;b = (1/ws)*ones(1,ws);a = 1;
%d_sst10 = filter(b,a,d_sst(:,ib2)')';
%d_pdsi10 = filter(b,a,d_pdsi(:,id2)')';
%x_sst10 = filter(b,a,x_sst(:,ir1)')';
%x_pdsi10 = filter(b,a,Xpdsi(:,ir1)')';

d_sst10 = d_sst(:,ib2);
d_pdsi10 = d_pdsi(:,id2);
x_sst10 = x_sst(:,ir1);
x_pdsi10 = Xpdsi(:,ir1);


% Compute the differences
dsz=size(d_sst10);
dszx=size(x_sst10);
df_sst10=zeros(dsz(1),dsz(2)-20);
df_pdsi10=zeros(dsz(1),dsz(2)-20);
xf_sst10=zeros(dszx(1),dszx(2)-20);
xf_pdsi10=zeros(dszx(1),dszx(2)-20);
dts1=zeros(4,length(id2)-20);
dts2=zeros(4,length(ir1)-20);
for i=1:(length(id2)-20)
   df_sst10(:,i)=mean(d_sst10(:,(i+10):(i+19)),2)-mean(d_sst10(:,i:(i+9)),2);
   df_pdsi10(:,i)=mean(d_pdsi10(:,(i+10):(i+19)),2)-mean(d_pdsi10(:,i:(i+9)),2);
   %df_sst10(:,i)=d_sst10(:,i+10)-d_sst10(:,i);
   %df_pdsi10(:,i)=d_pdsi10(:,i+10)-d_pdsi10(:,i);
   dts1(:,i)=[i+10,i+19,i,i+9]+analyrs2(1)-1;
end
for i=1:(length(ir1)-20)
   xf_sst10(:,i)=mean(x_sst10(:,(i+10):(i+19)),2)-mean(x_sst10(:,i:(i+9)),2);
   xf_pdsi10(:,i)=mean(x_pdsi10(:,(i+10):(i+19)),2)-mean(x_pdsi10(:,i:(i+9)),2);
   %xf_sst10(:,i)=x_sst10(:,i+10)-x_sst10(:,i);
   %xf_pdsi10(:,i)=x_pdsi10(:,i+10)-x_pdsi10(:,i);
   dts2(:,i)=[i+10,i+19,i,i+9]+analyrs1(1)-1;
end

% Find the largest decadal differences over the past few decades
Wd=W;Wd([1:10 86:96],:)=0;% zero out high latitudes
dfw=bsxfun(@times,df_sst10,Wd(:));
%dfw=bsxfun(@times,detrend(df_sst10')',Wd(:));
[dfmg,dfi]=sort(nansum(dfw,1));


load('./colormaps/MPL_BrBG.mat')
cmap1=cmap;
load('./colormaps/NCV_blue_red.mat')
s=load('coast');mlat=s.lat;mlon=s.long;


% Plot the largest differences, both warm and cold
figure
hold on
axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
i=1;varZ=reshape(df_sst10(:,dfi(i)),length(lat),length(lon));
pcolorm(lat,lon,varZ); % add additional value to avoid line in pcolor plots
plotm(mlat,mlon,'color',[0.5 0.5 0.5])
caxis([-0.75 0.75])
colormap(cmap) % colormap from 'WhiteYellowOrangeRed'
ttl=['HadISST ' num2str(dfi(i)+analyrs2(1)-1+10) '-' num2str(dfi(i)+analyrs2(1)-1+19) ' minus ' ...
num2str(dfi(i)+analyrs2(1)-1) '-' num2str(dfi(i)+analyrs2(1)-1+9) ];
title(ttl)
tightmap
ttl=strrep(ttl,' ','_'); % remove space
if prnt_fgs=='y';print('-dpng',['./figs/decdiff_' ttl '.png',]);end

figure
hold on
axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
i=1;varZ=reshape(df_pdsi10(:,dfi(i)),length(lat),length(lon));
pcolorm(lat,lon,varZ); % add additional value to avoid line in pcolor plots
plotm(mlat,mlon,'color',[0.5 0.5 0.5])
caxis([-1.25 1.25])
colormap(cmap1) % colormap from 'WhiteYellowOrangeRed'
ttl=['Dai PDSI ' num2str(dfi(i)+analyrs2(1)-1+10) '-' num2str(dfi(i)+analyrs2(1)-1+19) ' minus ' ...
num2str(dfi(i)+analyrs2(1)-1) '-' num2str(dfi(i)+analyrs2(1)-1+9) ];
title(ttl)
tightmap
ttl=strrep(ttl,' ','_'); % remove space
if prnt_fgs=='y';print('-dpng',['./figs/decdiff_' ttl '.png',]);end

figure
hold on
axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
i=length(dfi);varZ=reshape(df_sst10(:,dfi(i)),length(lat),length(lon));
pcolorm(lat,lon,varZ); % add additional value to avoid line in pcolor plots
plotm(mlat,mlon,'color',[0.5 0.5 0.5])
caxis([-0.75 0.75])
colormap(cmap) % colormap from 'WhiteYellowOrangeRed'
ttl=['HadISST ' num2str(dfi(i)+analyrs2(1)-1+10) '-' num2str(dfi(i)+analyrs2(1)-1+19) ' minus ' ...
num2str(dfi(i)+analyrs2(1)-1) '-' num2str(dfi(i)+analyrs2(1)-1+9) ];
title(ttl)
tightmap
ttl=strrep(ttl,' ','_'); % remove space
if prnt_fgs=='y';print('-dpng',['./figs/decdiff_' ttl '.png',]);end

figure
hold on
axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
i=length(dfi);varZ=reshape(df_pdsi10(:,dfi(i)),length(lat),length(lon));
pcolorm(lat,lon,varZ); % add additional value to avoid line in pcolor plots
plotm(mlat,mlon,'color',[0.5 0.5 0.5])
caxis([-1.25 1.25])
colormap(cmap1) % colormap from 'WhiteYellowOrangeRed'
ttl=['Dai PDSI ' num2str(dfi(i)+analyrs2(1)-1+10) '-' num2str(dfi(i)+analyrs2(1)-1+19) ' minus ' ...
num2str(dfi(i)+analyrs2(1)-1) '-' num2str(dfi(i)+analyrs2(1)-1+9)];
title(ttl)
tightmap
ttl=strrep(ttl,' ','_'); % remove space
if prnt_fgs=='y';print('-dpng',['./figs/decdiff_' ttl '.png',]);end



% Find the pattern correlation maps that are most similar
pcorr1=zeros(size(xf_sst10,2),1); 
pcorr2=zeros(size(xf_sst10,2),1); 
for j=1:size(xf_sst10,2)
   pcorr1(j)=nancorr(df_sst10(:,dfi(1)),xf_sst10(:,j));
   pcorr2(j)=nancorr(df_sst10(:,dfi(end)),xf_sst10(:,j));
end

% which are the highest correlations?
[pc1,ip1]=sort(pcorr1);
[pc2,ip2]=sort(pcorr2);


% Plot PHYDA with highest pattern correlations
figure
hold on
axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
i=length(ip1);varZ=reshape(xf_sst10(:,ip1(i)),length(lat),length(lon));
pcolorm(lat,lon,varZ); % add additional value to avoid line in pcolor plots
plotm(mlat,mlon,'color',[0.5 0.5 0.5])
caxis([-0.75 0.75])
colormap(cmap) % colormap from 'WhiteYellowOrangeRed'
ttl=['PHYDA SST ' num2str(ip1(i)+analyrs1(1)-1+10) '-' num2str(ip1(i)+analyrs1(1)-1+19) ' minus ' ...
num2str(ip1(i)+analyrs1(1)-1) '-' num2str(ip1(i)+analyrs1(1)-1+9)];
title(ttl)
tightmap
ttl=strrep(ttl,' ','_'); % remove space
if prnt_fgs=='y';print('-dpng',['./figs/decdiff_' ttl '.png',]);end

figure
hold on
axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
i=length(ip1);varZ=reshape(xf_pdsi10(:,ip1(i)),length(lat),length(lon));
pcolorm(lat,lon,varZ); % add additional value to avoid line in pcolor plots
plotm(mlat,mlon,'color',[0.5 0.5 0.5])
caxis([-1.25 1.25])
colormap(cmap1) % colormap from 'WhiteYellowOrangeRed'
ttl=['PHYDA PDSI ' num2str(ip1(i)+analyrs1(1)-1+10) '-' num2str(ip1(i)+analyrs1(1)-1+19) ' minus ' ...
num2str(ip1(i)+analyrs1(1)-1) '-' num2str(ip1(i)+analyrs1(1)-1+9)];
title(ttl)
tightmap
ttl=strrep(ttl,' ','_'); % remove space
if prnt_fgs=='y';print('-dpng',['./figs/decdiff_' ttl '.png',]);end


figure
hold on
axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
i=length(ip2);varZ=reshape(xf_sst10(:,ip2(i)),length(lat),length(lon));
pcolorm(lat,lon,varZ); % add additional value to avoid line in pcolor plots
plotm(mlat,mlon,'color',[0.5 0.5 0.5])
caxis([-0.75 0.75])
colormap(cmap) % colormap from 'WhiteYellowOrangeRed'
ttl=['PHYDA SST ' num2str(ip2(i)+analyrs1(1)-1+10) '-' num2str(ip2(i)+analyrs1(1)-1+19) ' minus ' ...
num2str(ip2(i)+analyrs1(1)-1) '-' num2str(ip2(i)+analyrs1(1)-1+9)];
title(ttl)
tightmap
ttl=strrep(ttl,' ','_'); % remove space
if prnt_fgs=='y';print('-dpng',['./figs/decdiff_' ttl '.png',]);end

figure
hold on
axesm('MapProjection','eqdcylin','MapLatLimit',[(min(lat)) (max(lat))],'MapLonLimit',[min(lon) max(lon)]);
i=length(ip2);varZ=reshape(xf_pdsi10(:,ip2(i)),length(lat),length(lon));
pcolorm(lat,lon,varZ); % add additional value to avoid line in pcolor plots
plotm(mlat,mlon,'color',[0.5 0.5 0.5])
caxis([-1.25 1.25])
colormap(cmap1) % colormap from 'WhiteYellowOrangeRed'
ttl=['PHYDA PDSI ' num2str(ip2(i)+analyrs1(1)-1+10) '-' num2str(ip2(i)+analyrs1(1)-1+19) ' minus ' ...
num2str(ip2(i)+analyrs1(1)-1) '-' num2str(ip2(i)+analyrs1(1)-1+9) ];
title(ttl)
tightmap
ttl=strrep(ttl,' ','_'); % remove space
if prnt_fgs=='y';print('-dpng',['./figs/decdiff_' ttl '.png',]);end































