
%% MAKE FIGURES COMPARING PDSI AND SPEI

varpth='/d1/nsteiger/climate-data/cesm-lme/';

load([varpth 'cesm_lme_009_pdsi_output_AWC_c_7_u2_gcm_09-Jun-2016_10:45:14.mat'],'pdsi_lndonly')
%load([varpth 'cesm_lme_009_spei_output_scl_15_krnl_e_u2_gcm_09-Jun-2016_14:42:00.mat'],'spei_lndonly','lndidx','lat','lon')
load([varpth 'cesm_lme_009_spei_output_scl_12_krnl_e_u2_gcm_09-Jun-2016_10:30:26.mat'],'spei_lndonly','lndidx','lat','lon')
%load([varpth 'cesm_lme_009_spei_output_scl_9_krnl_e_u2_gcm_09-Jun-2016_14:41:05.mat'],'spei_lndonly','lndidx','lat','lon')

% jja pdsi
pdsi_jja=mon2ann(pdsi_lndonly,6,8);
spei_jja=mon2ann(spei_lndonly,6,8);

% PLOT CORRELATIONS

d_corr=zeros(length(lndidx),1);
for i=1:length(lndidx)
    d_corr(i)=nancorr(pdsi_jja(i,:)',spei_jja(i,:)');
end

corr_f=NaN(length(lat)*length(lon),1);
corr_f(lndidx,:)=d_corr;
corr_f=reshape(corr_f,length(lat),length(lon));

%corr_f_1=corr_f;

%%

figure
hold on
%s=load('coast');mlat=s.lat; mlon=s.long;
axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
    'MapLonLimit',[min(lon) max(lon)]);
title('corr(PDSI,SPEI12), Gauss kernel')
pcolorm(lat,lon,corr_f)
colorbar 
caxis([0 1])
print('pdsi_spei_g_corr.eps','-depsc')

%%
figure
hold on
%s=load('coast');mlat=s.lat; mlon=s.long;
axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
    'MapLonLimit',[min(lon) max(lon)]);
title('corr(PDSI,SPEI12), Exp minus Gauss')
pcolorm(lat,lon,corr_f-corr_f_1)
colorbar 
%caxis([0 1])
%print('pdsi_spei_comparison_corr.png','-dpng')

%%
diff=corr_f-corr_f_1;

figure
hist(diff(:),100)
%print('pdsi_spei_comparison_corr_hist.png','-dpng')






%%

varpth='/d1/nsteiger/climate-data/gfdl-esm2m/';

load([varpth 'pdsi_output_AWC_7_u2_gcm_14-Apr-2016_14:59:06.mat'])
load([varpth 'spei_output_scl_18_krnl_g_u2_gcm_14-Apr-2016_14:53:16.mat'])

%%

figure
for i=36:72
hold on
s=load('coast');mlat=s.lat; mlon=s.long;
axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
    'MapLonLimit',[min(lon) max(lon)]);
title(['PDSI, month = ' num2str(i)])
pcolorm(lat,lon,pdsi_f(:,:,i))
colorbar 
caxis([-6 6])  
pause
end

%%

figure
for i=36:72
hold on
s=load('coast');mlat=s.lat; mlon=s.long;
axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
    'MapLonLimit',[min(lon) max(lon)]);
title(['SPEI, month = ' num2str(i)])
pcolorm(lat,lon,spei_f(:,:,i))
colorbar 
%caxis([-6 6])
pause
end

%% INDIVIDUAL TIME SERIES COMPARISONS
varpth='/d1/nsteiger/climate-data/gfdl-esm2m/';

load([varpth 'pdsi_output_AWC_7_u2_gcm_14-Apr-2016_14:59:06.mat'])

load([varpth 'spei_output_scl_18_krnl_g_u2_gcm_14-Apr-2016_14:53:16.mat'])
spei_lndonly_18=spei_lndonly;
load([varpth 'spei_output_scl_12_krnl_g_u2_gcm_14-Apr-2016_14:50:22.mat'])
spei_lndonly_12=spei_lndonly;
load([varpth 'spei_output_scl_9_krnl_g_u2_gcm_14-Apr-2016_14:49:17.mat'])
spei_lndonly_9=spei_lndonly;
load([varpth 'spei_output_scl_6_krnl_g_u2_gcm_14-Apr-2016_14:48:27.mat'])
spei_lndonly_6=spei_lndonly;
load([varpth 'spei_output_scl_3_krnl_g_u2_gcm_14-Apr-2016_14:48:11.mat'])
spei_lndonly_3=spei_lndonly;

%% INDIVIDUAL SITE EXAMPLE
j=100;

figure
hold on
plot(pdsi_lndonly(j,:))
plot(spei_lndonly_3(j,:))
plot(spei_lndonly_6(j,:))
plot(spei_lndonly_9(j,:))
plot(spei_lndonly_12(j,:))
plot(spei_lndonly_18(j,:))
legend('PDSI','SPEI_{3}','SPEI_{6}','SPEI_{9}','SPEI_{12}','SPEI_{18}')
title('Drought indices')
xlabel('Months')

%% INDVIDUAL SITE P AND PET

prlnd=pr_inmon(:,lndidx);
%%
j=478; 
figure
%for j=100:500
hold on
plot(pdsi_lndonly(j,:))
plot(prlnd(:,j))
legend('PDSI','PRECIP')
title(['PDSI & Precip ' num2str(j)])
box on
hold off
%pause
%end
%print('pdsi_full_dist.eps','-depsc')

%% PDSI DIST.

pdsiv=pdsi_lndonly(:);

figure
hist(pdsiv,400)
title('All PDSI')
print('pdsi_full_dist.eps','-depsc')

%% PDSI OUTLIERS

pdsiv=pdsi_lndonly(:);
out_PDSI=find(abs(pdsiv)>10);

figure
hist(pdsiv(out_PDSI),200)
title('Outlier PDSI')
print('pdsi_outlier_dist.eps','-depsc')

%% SPEI DIST.

speiv=spei_lndonly_18(:);

figure
hist(speiv,300)
title('All SPEI18')


%% PLOT CORRELATIONS

d_corr=zeros(length(lndidx),1);
for i=1:length(lndidx)
    d_corr(i)=nancorr(pdsi_lndonly(i,:)',spei_lndonly_12r(i,:)');
end

corr_f=NaN(length(lat)*length(lon),1);
corr_f(lndidx,:)=d_corr;
corr_f=reshape(corr_f,length(lat),length(lon));

%%

figure
hold on
%s=load('coast');mlat=s.lat; mlon=s.long;
axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
    'MapLonLimit',[min(lon) max(lon)]);
title('corr(PDSI,SPEI12), Rectangular kernel')
pcolorm(lat,lon,corr_f)
colorbar 
caxis([0.2 1])
print('pdsi_spei_r_corr.eps','-depsc')

%% CHECK DIFFERENCES WITH KERNELS

load('spei_output_scl_18_u2_gcm_14-Apr-2016_13:00:49.mat')
spei_lndonly_krnlg=spei_lndonly;
load('spei_output_scl_18_u2_gcm_14-Apr-2016_10:28:17.mat')
spei_lndonly_krnlr=spei_lndonly;

%%

for i=1:10
    figure
    hold on
   plot(spei_lndonly_krnlg(i,:))
   plot(spei_lndonly_krnlr(i,:))
    hold off
    legend('G','R')
end


%% PDSI FORMULATIONS COMPARISONS 

% NEED TO ADD ANALYSIS FOR MSE, BIAS, EVAR

clear
varpth='/d1/nsteiger/climate-data/gfdl-esm2m/';
load([varpth 'pdsi_output_AWC_7_u2_gcm_14-Apr-2016_14:59:06.mat'])
pdsi_bsln=pdsi_lndonly;
%load([varpth 'pdsi_output_AWC_7_u2_1ms_07-Apr-2016_14:56:18.mat'])
%pdsi_x=pdsi_lndonly;
%load([varpth 'pdsi_output_AWC_5_u2_gcm_07-Apr-2016_14:52:03.mat'])
%pdsi_x=pdsi_lndonly;
%load([varpth 'pdsi_output_AWC_5_u2_1ms_07-Apr-2016_11:10:20.mat'])
%pdsi_x=pdsi_lndonly;
%load([varpth 'pdsi_output_AWC_7_u10_gcm_15-Apr-2016_15:47:34.mat'])
%pdsi_x=pdsi_lndonly;
load([varpth 'pdsi_output_AWC_vrbl_u2_gcm_18-Apr-2016_13:17:01.mat'])
pdsi_x=pdsi_lndonly;

%%
% d_corr=zeros(length(lndidx),1);
% for i=1:length(lndidx)
%     d_corr(i)=nancorr(pdsi_bsln(i,:)',pdsi_x(i,:)');
% end
% 
% corr_f=NaN(length(lat)*length(lon),1);
% corr_f(lndidx,:)=d_corr;
% corr_f=reshape(corr_f,length(lat),length(lon));

d_bias=zeros(length(lndidx),1);
d_evar=zeros(length(lndidx),1);
d_mse=zeros(length(lndidx),1);
for i=1:length(lndidx)
    [d_mse(i),d_evar(i),d_bias(i)]=mse(pdsi_bsln(i,:),pdsi_x(i,:));
end

mse_f=NaN(length(lat)*length(lon),1);
mse_f(lndidx,:)=d_mse;
mse_f=reshape(mse_f,length(lat),length(lon));

evar_f=NaN(length(lat)*length(lon),1);
evar_f(lndidx,:)=d_evar;
evar_f=reshape(evar_f,length(lat),length(lon));

bias_f=NaN(length(lat)*length(lon),1);
bias_f(lndidx,:)=d_bias;
bias_f=reshape(bias_f,length(lat),length(lon));



%%
% figure
% hold on
% axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
%     'MapLonLimit',[min(lon) max(lon)]);
% title('corr PDSI(AWC7U2GCM), PDSI(AWC7U10)')
% pcolorm(lat,lon,corr_f)
% colorbar 
% caxis([0.5 1])
% print('corr_PDSI_AWC7U2GCM_PDSI_AWC7U10.eps','-depsc')

figure
hold on
axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
    'MapLonLimit',[min(lon) max(lon)]);
title('mse PDSI(AWC7 U2GCM), PDSI(AWCVar U2GCM)')
pcolorm(lat,lon,mse_f)
colorbar 
%caxis([0.5 1])
print('mse_PDSI_AWC7U2GCM_PDSI_AWCVarU2GCM.eps','-depsc')

figure
hold on
axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
    'MapLonLimit',[min(lon) max(lon)]);
title('evar PDSI(AWC7 U2GCM), PDSI(AWCVar U2GCM)')
pcolorm(lat,lon,evar_f)
colorbar 
%caxis([0.5 1])
print('evar_PDSI_AWC7U2GCM_PDSI_AWCVarU2GCM.eps','-depsc')

figure
hold on
axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
    'MapLonLimit',[min(lon) max(lon)]);
title('bias PDSI(AWC7 U2GCM), PDSI(AWCVar U2GCM)')
pcolorm(lat,lon,bias_f)
colorbar 
%caxis([0.5 1])
print('bias_PDSI_AWC7U2GCM_PDSI_AWCVarU2GCM.eps','-depsc')

%%

% Make PDSI figs

load('pdsi_output_05-Apr-2016_11:23:26.mat')
%%
% Reshape PDSI
pdsi_rshp=reshape(pdsi_locs,lyrs*12,length(lndidx))';
pdsi_f=NaN(length(lat)*length(lon),lyrs*12);
pdsi_f(lndidx,:)=pdsi_rshp;
pdsi_f=reshape(pdsi_f,length(lat),length(lon),lyrs*12);
pdsi_f=mon2ann(pdsi_f,1,12);
pdsi_f(pdsi_f==0)=NaN;
%%
varpth='/d1/nsteiger/climate-data/gfdl-esm2m/';
lat=ncread([varpth 'tas_Amon_GFDL-ESM2M_historical_r1i1p1_186101-200512.nc'],'lat');
lon=ncread([varpth 'tas_Amon_GFDL-ESM2M_historical_r1i1p1_186101-200512.nc'],'lon');
%%
figure
for i=1:145
hold on
s=load('coast');mlat=s.lat; mlon=s.long;
axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
    'MapLonLimit',[min(lon) max(lon)]);
title(['NOAA PDSI, month = ' num2str(i)])
pcolorm(lat,lon,jjapdsi_n(:,:,i))
colorbar 
caxis([-6 6])
pause
end
%%

% Figs for Sloan's data
figure
for i=1:145
hold on
s=load('coast');mlat=s.lat; mlon=s.long;
axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
    'MapLonLimit',[min(lon) max(lon)]);
title(['Sload PDSI, month = ' num2str(i)])
pcolorm(lat,lon,jjapdsi(:,:,i))
colorbar 
caxis([-6 6])
pause
end


%% COMPARE SLOAN'S CODE WITH MY NOAA CODE

i=138;

% jjapdsi_n=mon2ann(pdsi_f,6,8); %JJA PDSI 
% jjapdsi_n(jjapdsi_n==0)=NaN;

figure
hold on
s=load('coast');mlat=s.lat; mlon=s.long;
axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
    'MapLonLimit',[min(lon) max(lon)]);
title(['NOAA PDSI, year = ' num2str(1861+i-1)])
pcolorm(lat,lon,jjapdsi_n(:,:,i))
colorbar 
caxis([-6 6])

print(['pdsi_map_jjapdsi_noaa_yr_' num2str(1861+i-1) '.eps'],'-depsc')


% % land area
%varpth='/d1/nsteiger/climate-data/gfdl-esm2m/';
%sftlf=permute(ncread([varpth 'sftlf_fx_GFDL-ESM2M_historical_r0i0p0.nc'],'sftlf'),[2 1]);
%oceanidx=find((sftlf(:)<50));
% 
% jjapdsi=mil0001.jjapdsi;
% jjapdsi=reshape(permute(jjapdsi,[2 3 1]),length(lat)*length(lon),lyrs);
% 
% jjapdsi(oceanidx,:)=NaN;
% 
% jjapdsi=reshape(jjapdsi,length(lat),length(lon),lyrs);


figure
hold on
s=load('coast');mlat=s.lat; mlon=s.long;
axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
    'MapLonLimit',[min(lon) max(lon)]);
title(['Sloan PDSI, year = ' num2str(1861+i-1)])
pcolorm(lat,lon,jjapdsi(:,:,i))
colorbar 
caxis([-6 6])

print(['pdsi_map_jjapdsi_sloan_yr_' num2str(1861+i-1) '.eps'],'-depsc')

%% FIND WHERE OUTLIERS ACTUALLY ARE

out_noaa=find(abs(jjapdsi_n(:))>10);
out_sloan=find(abs(jjapdsi(:))>10);

%%

pdsi_jja=mon2ann(pdsi_f,6,8);
pdsi_jja(pdsi_jja==0)=NaN;
howmany=zeros(lyrs,1);
for j=1:lyrs
   [rw_o,cl_o]=find(abs(pdsi_jja(:,:,j))>10); 
   howmany(j)=length(rw_o);
end

%%

j=10;
figure
hold on
s=load('coast');mlat=s.lat; mlon=s.long;
axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
    'MapLonLimit',[min(lon) max(lon)]);
title(['NOAA PDSI, yr = ' num2str(1861+j-1)],'Interpreter','none')
pcolorm(lat,lon,pdsi_jja(:,:,j))
colorbar 
print(['pdsi_map_jjapdsi_sloan_yr_' num2str(1861+j-1) '.eps'],'-depsc')

%% FIND OUT HOW PARAMETER CHOICES IN PDSI CALCULATION CHANGE RESULTS
clear
varpth='/d1/nsteiger/climate-data/gfdl-esm2m/';

%load([varpth 'pdsi_output_AWC_5_u2_1ms_07-Apr-2016_11:10:20.mat'])
load([varpth 'pdsi_output_AWC_5_u2_gcm_07-Apr-2016_14:52:03.mat'])
%load([varpth 'pdsi_output_AWC_7_u2_1ms_07-Apr-2016_14:56:18.mat'])
%load([varpth 'pdsi_output_AWC_7_u2_gcm_05-Apr-2016_16:08:22.mat'])
pdsi_jja_1=mon2ann(pdsi_f,6,8); %JJA PDSI 
pdsi_jja_1(pdsi_jja_1==0)=NaN;


%load([varpth 'pdsi_output_AWC_5_u2_1ms_07-Apr-2016_11:10:20.mat'])
%load([varpth 'pdsi_output_AWC_5_u2_gcm_07-Apr-2016_14:52:03.mat'])
load([varpth 'pdsi_output_AWC_7_u2_1ms_07-Apr-2016_14:56:18.mat'])
%load([varpth 'pdsi_output_AWC_7_u2_gcm_05-Apr-2016_16:08:22.mat'])
pdsi_jja_2=mon2ann(pdsi_f,6,8); %JJA PDSI 
pdsi_jja_2(pdsi_jja_2==0)=NaN;


% mean squared error difference
[MSE,EVAR,BIAS]= mse(pdsi_jja_1,pdsi_jja_2);


% mean state differences
m1=nanmean(pdsi_jja_1,3);
m1(m1==0)=NaN;
m2=nanmean(pdsi_jja_2,3);
m2(m2==0)=NaN;


figure
hold on
s=load('coast');mlat=s.lat; mlon=s.long;
axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
    'MapLonLimit',[min(lon) max(lon)]);
title(['NOAA PDSI AWC_5_u2_1ms'],'Interpreter','none')
pcolorm(lat,lon,m1)
colorbar 
%caxis([-6 6])

figure
hold on
s=load('coast');mlat=s.lat; mlon=s.long;
axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
    'MapLonLimit',[min(lon) max(lon)]);
title(['Sloan PDSI AWC_7_u2_gcm'],'Interpreter','none')
pcolorm(lat,lon,m2)
colorbar 
%caxis([-6 6])


%% INDIVIDUAL YEAR DIFFERENCES
i=50;
figure
hold on
s=load('coast');mlat=s.lat; mlon=s.long;
axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
    'MapLonLimit',[min(lon) max(lon)]);
title(['NOAA PDSI AWC_5_u2_gcm'],'Interpreter','none')
pcolorm(lat,lon,pdsi_jja_1(:,:,i))
colorbar 
%caxis([-6 6])

figure
hold on
s=load('coast');mlat=s.lat; mlon=s.long;
axesm('MapProjection','eqdcylin','MapLatLimit',[min(lat) max(lat)],...
    'MapLonLimit',[min(lon) max(lon)]);
title(['Sloan PDSI AWC_7_u2_1ms'],'Interpreter','none')
pcolorm(lat,lon,pdsi_jja_2(:,:,i))
colorbar 
%caxis([-6 6])

%% PDSI OUTLIERS: WHY?



% CHECK SOIL MOISTURE
varpth='/d1/nsteiger/climate-data/gfdl-esm2m/';
syr=1861;
eyr=2005;
lyrs=length(syr:eyr);

load([varpth 'pdsi_output_AWC_7_u2_gcm_11-Apr-2016_15:36:57.mat']);

lat=ncread([varpth 'tas_Amon_GFDL-ESM2M_historical_r1i1p1_186101-200512.nc'],'lat');
lon=ncread([varpth 'tas_Amon_GFDL-ESM2M_historical_r1i1p1_186101-200512.nc'],'lon');

mrso=permute(ncread([varpth 'mrso_Lmon_GFDL-ESM2M_historical_r1i1p1_186101-200512.nc'],'mrso'),[2 1 3]);
pr=permute(ncread([varpth 'pr_Amon_GFDL-ESM2M_historical_r1i1p1_186101-200512.nc'],'pr'),[2 1 3]);
sftlf=permute(ncread([varpth 'sftlf_fx_GFDL-ESM2M_historical_r0i0p0.nc'],'sftlf'),[2 1]);

lndidx=find((sftlf>50)); % get grid cells with < 50% land
ocnidx=find((sftlf<50)); 

% find outliers
pdsi_jja=mon2ann(pdsi_f,6,8);
pdsi_jja(pdsi_jja==0)=NaN;
howmany=zeros(lyrs,1);
rwo=[];clo=[];
for j=1:lyrs
   [rw_o,cl_o]=find(abs(pdsi_jja(:,:,j))>10); 
   howmany(j)=length(rw_o);
   rwo(1:length(rw_o),j)=rw_o;
   clo(1:length(cl_o),j)=cl_o;
end

mrso_jja=mon2ann(mrso,6,8);
mrso_jja(mrso_jja==0)=NaN;

pr_jja=mon2ann(pr,6,8);
pr_jja(pr_jja==0)=NaN;

pr_jjar=reshape(pr_jja,length(lat)*length(lon),lyrs);
pr_jjar(ocnidx,:)=NaN;
pr_jja=reshape(pr_jjar,length(lat),length(lon),lyrs);

%%
close all

i=36;
figure
pcolor(mrso_jja(:,:,i))
title('mrso')
%%
grtzr=find(howmany>0)';
for i=grtzr(1:20)
figure
pcolor(pr_jja(:,:,i))
title(['pr ' num2str(i)])
colorbar
end
%%
grtzr=find(howmany>0)';
for i=grtzr(end-10:end)
figure
pcolor(pdsi_jja(:,:,i))
title(['pdsi ' num2str(i)])
colorbar
end

%%

% pdsi 142, row 56, col 11












