%%

clc
clear 
close 

addpath('/Users/yang560/Dropbox/PNNL/Matlab_functions/utilities/')


%% 2     3       4     5        6         7      8     9     10     11
% 'BB  DRYSMC   F11  MAXSMC   REFSMC   SATPSI  SATDK SATDW WLTSMC  QTZ;
sct = ncread('geo_em.d04.nc.gssurgo','SCT_DOM');

lat2d = ncread('geo_em.d04.nc.gssurgo', 'XLAT_M');
lon2d = ncread('geo_em.d04.nc.gssurgo', 'XLONG_M');
hgt = ncread('geo_em.d04.nc.gssurgo', 'HGT_M');

SMPara_old = load('SoilTBL_OLD.txt'); % 2nd BB; 5th maxsmc; 8th SATDK;
SMPara_new = load('SoilTBL_NEW.txt'); % 2nd BB; 5th maxsmc; 8th SATDK;

SMPara = cat(3, SMPara_old, SMPara_new);

for ii = 1:392
    for jj = 1:392
        bb(:,ii,jj)      = SMPara(sct(ii,jj),2,:);
        drysmc(:,ii,jj)  = SMPara(sct(ii,jj),3,:);
        f11(:,ii, jj)    = SMPara(sct(ii,jj),4,:);
        maxsmc(:,ii, jj) = SMPara(sct(ii,jj),5,:);
        refsmc(:,ii, jj) = SMPara(sct(ii,jj),6,:);
        satpsi(:,ii, jj) = SMPara(sct(ii,jj),7,:);
        satdk(:,ii,jj)   = SMPara(sct(ii,jj),8,:);
        satdw(:,ii,jj)   = SMPara(sct(ii,jj),9,:);
        wltsmc(:,ii,jj)  = SMPara(sct(ii,jj),10,:);
        qtz(:,ii,jj)     = SMPara(sct(ii,jj),11,:);
    end 
end

bb_diff      = squeeze(bb(2,:,:)     - bb(1,:,:));
drysmc_diff  = squeeze(drysmc(2,:,:) - drysmc(1,:,:));
f11_diff     = squeeze(f11(2,:,:)    - f11(1,:,:));
maxsmc_diff  = squeeze(maxsmc(2,:,:) - maxsmc(1,:,:));
refsmc_diff  = squeeze(refsmc(2,:,:) - refsmc(1,:,:));
satpsi_diff  = squeeze(satpsi(2,:,:) - satpsi(1,:,:));
satdk_diff   = squeeze(satdk(2,:,:)  - satdk(1,:,:));
satdw_diff   = squeeze(satdw(2,:,:)  - satdw(1,:,:));
wltsmc_diff  = squeeze(wltsmc(2,:,:) - wltsmc(1,:,:));
qtz_diff     = squeeze(qtz(2,:,:)    - qtz(1,:,:));


%% plot_section 
close all
make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.04 0.01], [0.06 0.025], [0.03 0.03]);
if ~make_it_tight,  clear subplot;  end

set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [0 0 1000 600 ])
% figure default 
figure(1)
%---------------------------------------------------------------
subplot(3,3,1)
surf(lon2d, lat2d,bb_diff); view(2); shading interp;
caxis([ -1 1]);colormap (bluewhitered(20)); 
h=colorbar; title(h, '[-]')
title_left('a) Beta ')
ylim([35.1, 38.1])
xlim([-99.3,-95.6])
set(gca,'xticklabel',{[]})
box on;
%---------------------------------------------------------------
subplot(3,3,2)
surf(lon2d, lat2d,drysmc_diff); view(2); shading interp;
max_lim = max(abs(drysmc_diff(:)));
caxis([-1*max_lim max_lim]);
colormap (bluewhitered);
h=colorbar; title(h, '[m^3/m^{3}]')
title_left('b) Dry soil moisture thershold  ')
ylim([35.1, 38.1])
xlim([-99.3,-95.6])
set(gca,'xticklabel',{[]})
set(gca,'yticklabel',{[]})
box on;

%---------------------------------------------------------------
subplot(3,3,3)
surf(lon2d, lat2d,maxsmc_diff); view(2); shading interp;
max_lim = max(abs(maxsmc_diff(:)));
caxis([-1*max_lim max_lim]);
colormap (bluewhitered);
h=colorbar; title(h, '[m^3/m^{3}]')
title_left('c) Porosity  ')
ylim([35.1, 38.1])
xlim([-99.3,-95.6])
set(gca,'xticklabel',{[]})
set(gca,'yticklabel',{[]})
box on;

%---------------------------------------------------------------
subplot(3,3,4)
surf(lon2d, lat2d,refsmc_diff); view(2); shading interp;
max_lim = max(abs(refsmc_diff(:)));
caxis([-1*max_lim max_lim]);
colormap (bluewhitered);
h=colorbar; title(h, '[m^3/m^{3}]')
title_left('d) Field capacity ')
ylim([35.1, 38.1])
xlim([-99.3,-95.6])
set(gca,'xticklabel',{[]})
box on;

%---------------------------------------------------------------
subplot(3,3,5)
surf(lon2d, lat2d,satpsi_diff); view(2); shading interp;
max_lim = max(abs(satpsi_diff(:)));
caxis([-1*max_lim max_lim]);
colormap (bluewhitered);
h=colorbar; title(h, '[m]' )
title_left('e) Soil potential')
ylim([35.1, 38.1])
xlim([-99.3,-95.6])
set(gca,'xticklabel',{[]})
set(gca,'yticklabel',{[]})
box on;

%---------------------------------------------------------------
subplot(3,3,6)
surf(lon2d, lat2d,satdk_diff); view(2); shading interp;
max_lim = max(abs(satdk_diff(:)));
caxis([-1*max_lim max_lim]);
colormap (bluewhitered);
h=colorbar; title(h, '[m/s]')
title_left('f) Saturated hydraulic conductivity  ')
ylim([35.1, 38.1])
xlim([-99.3,-95.6])
set(gca,'yticklabel',{[]})
box on;

%---------------------------------------------------------------
subplot(3,3,7)
surf(lon2d, lat2d,satdw_diff); view(2); shading interp;
max_lim = max(abs(satdw_diff(:)));
caxis([-1*max_lim max_lim]);
colormap (bluewhitered);
h=colorbar; title(h, '[m^2/s]')
title_left('g) Saturated hydraulic diffusivity  ')
ylim([35.1, 38.1])
xlim([-99.3,-95.6])
box on;

%---------------------------------------------------------------
subplot(3,3,8)
surf(lon2d, lat2d,wltsmc_diff); view(2); shading interp;
max_lim = max(abs(wltsmc_diff(:)));
caxis([-1*max_lim max_lim]);
colormap (bluewhitered);
h=colorbar; title(h, '[m^3/m^{3}]')
title_left('h) Wilting point soil moisture  ')
ylim([35.1, 38.1])
xlim([-99.3,-95.6])
set(gca,'yticklabel',{[]})
box on;

%---------------------------------------------------------------

print('-depsc','SOIL_PARM_COMBINED')


%% PLOT SEASONAL SCATTER PLOT

% calculate different parameter sensitivity to SM1 and SM

diff_SM1 = ncread('../../16_soil_type_diff/DATA_DIFF/diff_SM1_1-9.nc', 'diff_SOIL_M1');
diff_SM2 = ncread('../../16_soil_type_diff/DATA_DIFF/diff_SM2_1-9.nc', 'diff_SOIL_M2');
diff_SM3 = ncread('../../16_soil_type_diff/DATA_DIFF/diff_SM3_1-9.nc', 'diff_SOIL_M3');
diff_SM4 = ncread('../../16_soil_type_diff/DATA_DIFF/diff_SM4_1-9.nc', 'diff_SOIL_M4');

diff_SM = diff_SM1*0.1 + diff_SM2*0.3 + diff_SM3*0.6 + diff_SM4*1.0;

diff_SM1_smtbl = squeeze(diff_SM1(:,:,2,:));   % 2 --- CORRESPONDS TO ROPT3-ROPT2
diff_SM2_smtbl = squeeze(diff_SM2(:,:,2,:));   % 2 --- CORRESPONDS TO ROPT3-ROPT2
diff_SM3_smtbl = squeeze(diff_SM3(:,:,2,:));   % 2 --- CORRESPONDS TO ROPT3-ROPT2
diff_SM4_smtbl = squeeze(diff_SM4(:,:,2,:));   % 2 --- CORRESPONDS TO ROPT3-ROPT2

diff_SM_smtbl  = squeeze(diff_SM(:,:,2,:)); 

% close all
% plot_scatter_xy(maxsmc_diff, diff_SM1_smtbl); suptitle('maxsmc'); print('-depsc', 'maxsmc_scatter_season') 
% plot_scatter_xy(satdk_diff,  diff_SM1_smtbl); suptitle('satdk'); print('-depsc', 'satdk_scatter_season') 
% plot_scatter_xy(drysmc_diff, diff_SM1_smtbl); suptitle('drysmc'); print('-depsc', 'drysmc_scatter_season') 
% plot_scatter_xy(bb_diff,     diff_SM1_smtbl);     suptitle('bb'); print('-depsc', 'bb_scatter_season') 
% plot_scatter_xy(refsmc_diff, diff_SM1_smtbl); suptitle('refsmc'); print('-depsc', 'refsmc_scatter_season') 
% plot_scatter_xy(wltsmc_diff, diff_SM1_smtbl); suptitle('wltsmc'); print('-depsc', 'wltsmc_scatter_season') 
% plot_scatter_xy(satdw_diff,  diff_SM1_smtbl);  suptitle('satdw'); print('-depsc', 'satdw_scatter_season') 
% plot_scatter_xy(satpsi_diff, diff_SM1_smtbl); suptitle('satpsi'); print('-depsc', 'satpsi_scatter_season') 

%% PLOT ANNUAL SCATTER PLOT....

close all

make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.04 0.02], [0.06 0.025], [0.03 0.03]);
if ~make_it_tight,  clear subplot;  end

figure
set(gcf,'PaperPositionMode', 'auto');
set(gcf,'Position', [0 0 1000 600])

subplot(3,3,1);plot_scatter_xy_annual(maxsmc_diff, diff_SM_smtbl,'SM', 'maxsmc'); title_left('a) maxsmc'); %print('-depsc', 'maxsmc_scatter_annual') 
subplot(3,3,2);plot_scatter_xy_annual(satdk_diff,  diff_SM_smtbl,'SM', 'satdk' ); title_left('b) satdk');  %print('-depsc', 'satdk_scatter_annual') 
subplot(3,3,3);plot_scatter_xy_annual(drysmc_diff, diff_SM_smtbl,'SM', 'drysmc'); title_left('c) drysmc'); %print('-depsc', 'drysmc_scatter_annual') 
subplot(3,3,4);plot_scatter_xy_annual(bb_diff,     diff_SM_smtbl,'SM', 'bb'    ); title_left('d) bb');     %print('-depsc', 'bb_scatter_annual')  
subplot(3,3,5);plot_scatter_xy_annual(refsmc_diff, diff_SM_smtbl,'SM', 'refsmc'); title_left('e) refsmc'); %print('-depsc', 'refsmc_scatter_annual') 
subplot(3,3,6);plot_scatter_xy_annual(wltsmc_diff, diff_SM_smtbl,'SM', 'wltsmc'); title_left('f) wltsmc'); %print('-depsc', 'wltsmc_scatter_annual') 
subplot(3,3,7);plot_scatter_xy_annual(satdw_diff,  diff_SM_smtbl,'SM', 'satdw' ); title_left('g) satdw');  %print('-depsc', 'satdw_scatter_annual') 
subplot(3,3,8);plot_scatter_xy_annual(satpsi_diff, diff_SM_smtbl,'SM', 'satpsi'); title_left('h) satpsi'); %print('-depsc', 'satpsi_scatter_annual') 
print('-depsc', 'SM_TBL_SCATTER_ANNUAL') 

%% PLOT spatial map for different soil category 1D
close all

make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.075 0.045], [0.06 0.04], [0.05 0.12]);
if ~make_it_tight,  clear subplot;  end

soil_type = {'sand','loamy sand','sandy loam','silt loam',...
    'silt','loam','sandy clay loam','silty clay loam',...
    'clay loam','sandy clay','silty clay','clay',...
    'organic material','water','bedrock','other'};
season = {'a) DJF', 'b) MAM', 'c) JJA', 'd) SON'};

for ii = 1:16
    
   figure
   set(gcf,'PaperPositionMode', 'auto');
    set(gcf,'Position', [0 0 600 400])

   for isn = 1:4
   diff_SM_single  = squeeze(diff_SM_smtbl(:,:,isn));
   diff_SM_single(sct~=ii) = nan;
   subplot(2,2,isn)
   plot_sct_diff_seasonal(lat2d, lon2d, diff_SM_single);
   title_left(season(isn))
   text (-97.5, 38.25, soil_type(ii))

   if isn== 4
   hp4 = get(subplot(2,2,4),'Position');
   h=colorbar('Position',[hp4(1)+hp4(3)+0.04  hp4(2)  0.02  hp4(2)+hp4(3)*1.95]);
   title(h, '\DeltaSM [m]')
   end
   
   end
   
   save_fig_name = [num2str(ii), '_', char(soil_type(ii)),'_seasonal'];
   print('-depsc', save_fig_name);
   
end

%% PLOT spatial map for different soil category 3D
close all

make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.045 0.075], [0.06 0.04], [0.05 0.09]);
if ~make_it_tight,  clear subplot;  end

soil_type = {'sand','loamy sand','sandy loam','silt loam',...
    'silt','loam','sandy clay loam','silty clay loam',...
    'clay loam','sandy clay','silty clay','clay',...
    'organic material','water','bedrock','other'};
season = {'a) DJF', 'b) MAM', 'c) JJA', 'd) SON'};

for ii = 1:16
    
   figure
   set(gcf,'PaperPositionMode', 'auto');
    set(gcf,'Position', [0 0 700 400])

   for isn = 1:4
   diff_SM_single  = squeeze(diff_SM_smtbl(:,:,isn));
   diff_SM_single(sct~=ii) = nan;
   subplot(2,2,isn)
   plot_sct_3d_diff_seasonal(lat2d, lon2d, hgt, diff_SM_single);
   title_left(season(isn))
   text (-97.5, 38.,600, soil_type(ii))

   if isn== 4
   hp4 = get(subplot(2,2,4),'Position');
   h=colorbar('Position', [hp4(1)+hp4(3)+0.03  hp4(2)  0.02  hp4(2)+hp4(3)*1.95]);
   title(h, '\DeltaSM [m]')
   end
   
   end
   
%    suptitle(soil_type(ii))

   save_fig_name = [num2str(ii), '_3D_', char(soil_type(ii)),'_seasonal'];
   print('-depsc', save_fig_name);
   
end

%% ---- PLOT ANNUAL SPATIAL PATTERN ----
close all
make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.04 0.01], [0.06 0.04], [0.03 0.07]);
if ~make_it_tight,  clear subplot;  end
ann = char(97:122);

figure
set(gcf,'PaperPositionMode', 'auto');
set(gcf,'Position', [0 0 1000 600])
for ii = 1:12
   
   subplot(3,4,ii)
   diff_SM1_single  = mean(diff_SM_smtbl,3);
   diff_SM1_single(sct~=ii) = nan;
   plot_sct_diff(lat2d, lon2d, diff_SM1_single); 
   title_left([ann(ii), ') ',char(soil_type(ii))])
   
   if (ii==1 || ii==2 ||ii==3 ||ii==4 ||ii==5 ||ii==6 ||ii==7 ||ii==8 )
       set(gca,'xticklabel',{[]})      
   end
   
   if ( ii==2 ||ii==3 ||ii==4 ||ii==6 ||ii==7 ||ii==8 ||ii==10 ||ii==11 ||ii==12 )
       set(gca,'yticklabel',{[]}) 
   end
   
   hp4 = get(subplot(3,4,12),'Position');

   h=colorbar('Position', [hp4(1)+hp4(3)+0.02  hp4(2)+0.08  0.01  hp4(2)+hp4(3)*3.1]);
   colormap (bluewhitered(20));
   title(h, '\DeltaSM [m]')
end
save_fig_name = ['SM_TBL_SM_DIFF_SPATIAL_PATTERN'];
print('-depsc', save_fig_name);

%% ---- PLOT BAR PLOT ----

[diff_stats_SM1] = func_calc_soil_bar(sct,diff_SM1_smtbl);
[diff_stats_SM2] = func_calc_soil_bar(sct,diff_SM2_smtbl);
[diff_stats_SM3] = func_calc_soil_bar(sct,diff_SM3_smtbl);
[diff_stats_SM4] = func_calc_soil_bar(sct,diff_SM4_smtbl);
[diff_stats_SM ] = func_calc_soil_bar(sct,diff_SM_smtbl);

%
close all
make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.01 0.01], [0.12 0.04], [0.12 0.02]);
if ~make_it_tight,  clear subplot;  end
ann = char(97:122);

figure
set(gcf,'PaperPositionMode', 'auto');
set(gcf,'Position', [0 0 400 600])
subplot(4,1,1)
plot_bar_with_errorbar(diff_stats_SM1);set(gca,'xticklabel',{[]});  text( 0.2, 0.06, 'a) SM1'); ylabel('SM diff. [m^3 m^{-3}]')
subplot(4,1,2)
plot_bar_with_errorbar(diff_stats_SM2);set(gca,'xticklabel',{[]})  ; text( 0.2, 0.06, 'b) SM2'); ylabel('SM diff. [m^3 m^{-3}]')
subplot(4,1,3)
plot_bar_with_errorbar(diff_stats_SM3);set(gca,'xticklabel',{[]})  ; text( 0.2, 0.06, 'c) SM3'); ylabel('SM diff. [m^3 m^{-3}]')
subplot(4,1,4)
plot_bar_with_errorbar(diff_stats_SM4);   text( 0.2, 0.06, 'd) SM4')  ;ylabel('SM diff. [m^3 m^{-3}]')
xticks([1:16]);
xticklabels(soil_type); 
xtickangle(30);

save_fig_name = 'SM_TBL_SM_DIFF_BAR_PLOT';
print('-depsc', save_fig_name);

%%
figure
set(gcf,'PaperPositionMode', 'auto');
set(gcf,'Position', [0 0 400 200])
plot_bar_with_errorbar(diff_stats_SM);set(gca,'xticklabel',{[]}); 
title_left( 'SM\_2m '); ylim([-0.1 0.1]);
ylabel('SM diff. [m]')
save diff_stats_SM.mat diff_stats_SM

save_fig_name = 'SM_TBL_SM_total_DIFF_BAR_PLOT';
print('-depsc', save_fig_name);

