%%

clc
clear 
close all
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');

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 SEASONAL SCATTER PLOT

% calculate different parameter sensitivity to LH and SM

diff_LH = ncread('../../16_soil_type_diff/DATA_DIFF/diff_LH_1-9.nc', 'diff_LH');
diff_SH = ncread('../../16_soil_type_diff/DATA_DIFF/diff_SH_1-9.nc', 'diff_SH');

diff_LH_smtbl = squeeze(diff_LH(:,:,2,:));   % 2 --- CORRESPONDS TO ROPT3-ROPT2
diff_SH_smtbl = squeeze(diff_SH(:,:,2,:));

% close all
% plot_scatter_xy(maxsmc_diff, diff_LH_smtbl); suptitle('maxsmc'); print('-depsc', 'maxsmc_scatter_season') 
% plot_scatter_xy(satdk_diff,  diff_LH_smtbl); suptitle('satdk'); print('-depsc', 'satdk_scatter_season') 
% plot_scatter_xy(drysmc_diff, diff_LH_smtbl); suptitle('drysmc'); print('-depsc', 'drysmc_scatter_season') 
% plot_scatter_xy(bb_diff, diff_LH_smtbl);     suptitle('bb'); print('-depsc', 'bb_scatter_season') 
% plot_scatter_xy(refsmc_diff, diff_LH_smtbl); suptitle('refsmc'); print('-depsc', 'refsmc_scatter_season') 
% plot_scatter_xy(wltsmc_diff, diff_LH_smtbl); suptitle('wltsmc'); print('-depsc', 'wltsmc_scatter_season') 
% plot_scatter_xy(satdw_diff, diff_LH_smtbl);  suptitle('satdw'); print('-depsc', 'satdw_scatter_season') 
% plot_scatter_xy(satpsi_diff, diff_LH_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_LH_smtbl,'LH', 'maxsmc'); title_left('a) maxsmc'); %print('-depsc', 'maxsmc_scatter_annual') 
subplot(3,3,2);plot_scatter_xy_annual(satdk_diff,  diff_LH_smtbl,'LH', 'satdk' ); title_left('b) satdk');  %print('-depsc', 'satdk_scatter_annual') 
subplot(3,3,3);plot_scatter_xy_annual(drysmc_diff, diff_LH_smtbl,'LH', 'drysmc'); title_left('c) drysmc'); %print('-depsc', 'drysmc_scatter_annual') 
subplot(3,3,4);plot_scatter_xy_annual(bb_diff,     diff_LH_smtbl,'LH', 'bb'    ); title_left('d) bb');     %print('-depsc', 'bb_scatter_annual')  
subplot(3,3,5);plot_scatter_xy_annual(refsmc_diff, diff_LH_smtbl,'LH', 'refsmc'); title_left('e) refsmc'); %print('-depsc', 'refsmc_scatter_annual') 
subplot(3,3,6);plot_scatter_xy_annual(wltsmc_diff, diff_LH_smtbl,'LH', 'wltsmc'); title_left('f) wltsmc'); %print('-depsc', 'wltsmc_scatter_annual') 
subplot(3,3,7);plot_scatter_xy_annual(satdw_diff,  diff_LH_smtbl,'LH', 'satdw' ); title_left('g) satdw');  %print('-depsc', 'satdw_scatter_annual') 
subplot(3,3,8);plot_scatter_xy_annual(satpsi_diff, diff_LH_smtbl,'LH', 'satpsi'); title_left('h) satpsi'); %print('-depsc', 'satpsi_scatter_annual') 

print('-depsc', 'SM_TBL_LH_SCATTER_ANNUAL') 
%% PLOT spatial map for different soil category
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_LH_single  = squeeze(diff_LH_smtbl(:,:,isn));
   diff_LH_single(sct~=ii) = nan;
   subplot(2,2,isn)
   plot_sct_diff_seasonal(lat2d, lon2d, diff_LH_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, '\DeltaLH[W m^{-2}]')
   end
   
   end
   
%    suptitle(soil_type(ii))

   save_fig_name = [num2str(ii), '_', 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_LH_single  = mean(diff_LH_smtbl,3);
   diff_LH_single(sct~=ii) = nan;
   plot_sct_diff_annual(lat2d, lon2d, diff_LH_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.025  hp4(2)+0.08  0.01  hp4(2)+hp4(3)*3.1])
    colormap (bluewhitered(20));
    title(h,'\DeltaLH [W m^{-2}]')
end

save_fig_name = ['SM_TBL_LH_DIFF_SPATIAL_PATTERN'];

print('-depsc', save_fig_name);

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

[diff_stats_LH] = func_calc_soil_bar(sct,diff_LH_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 200])

plot_bar_with_errorbar(diff_stats_LH);
ylim ([-6 10]); 
ylabel('LH diff. (W m^{-2})')
xticks([1:16]);
xticklabels(soil_type); 
xtickangle(45);

save diff_stats_LH.mat diff_stats_LH

% x=1:16;
% data    = diff_stats(:,1);
% errlow  = diff_stats(:,2)* -1.;
% errhigh = diff_stats(:,2);
% 
% % Find where data is positive or negative.
% positiveIndexes = data >= 0;
% negativeIndexes = data <  0;
% % Plot column plot (bar chart).
% y1 = data;
% y1(~positiveIndexes) = nan;
% y2 = data; 
% y2(~negativeIndexes) = nan;
% 
% bar(x, y1, 'r')
% hold on
% bar(x, y2, 'b')
% hold on
% 
% er = errorbar(x,data,errlow,errhigh);    
% er.Color = [0 0 0];                            
% er.LineStyle = 'none';  
% 
% hold off
% grid on;



print('-depsc', 'Categorical_LH_diff')

%%