%% look at the difference in soil type between gssurgo and default.

clc
clear 
close all 
addpath('/Users/yang560/Dropbox/PNNL/Matlab_functions/utilities/')

%% 
soilt_gssurgo = ncread('geo_em.d04.nc.gssurgo', 'SCT_DOM');
soilt_default = ncread('geo_em.d04.nc.backup', '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');

lat2d = double (lat2d);
lon2d = double (lon2d);
hgt   = double (hgt);

%%
figure(1)
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 600 300])
surf(lon2d, lat2d, hgt,soilt_gssurgo); 
shading interp; view(20,50); 
xlabel('Longitude')
ylabel('Latitude')
zlabel('Surface Elevation (m)')
colormap (jet(16))
caxis([.5, 16.5])
h = colorbar;
h.Ticks = linspace(1,16,16);
h.TickLabels = {'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'};

title(h,'Soil Type')
zlim([0 800]);
print('-depsc', 'soiltype_3D_gssurgo')

figure(2)
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 600 300])
surf(lon2d, lat2d, hgt,soilt_gssurgo); 
shading interp; view(2); 
xlabel('Longitude')
ylabel('Latitude')
zlabel('Surface Elevation (m)')
colormap (jet(16))
caxis([.5, 16.5])
h = colorbar
h.Ticks = linspace(1,16,16);
h.TickLabels = {'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'};
title(h,'Soil Type')
zlim([0 800]);
print('-depsc', 'soiltype_1D_gssurgo')

%%
figure(3)
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 600 300])
surf(lon2d, lat2d, hgt,soilt_default); 
shading interp; view(20,50); 
xlabel('Longitude')
ylabel('Latitude')
zlabel('Surface Elevation (m)')
colormap (jet(16))
caxis([.5,16.5]) % value from the exact bottom to the exact top of the colorbar
h = colorbar
h.Ticks = linspace(1,16,16);
h.TickLabels = {'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'};

title(h,'Soil Type')
zlim([0 800]);
print('-depsc', 'soiltype_3D_default')

%%
figure(4)
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 600 300])
surf(lon2d, lat2d, hgt,soilt_default); 
shading interp; view(2); 
xlabel('Longitude')
ylabel('Latitude')
zlabel('Surface Elevation (m)')
colormap (jet(16))
caxis([.5, 16.5])
h = colorbar;
h.Ticks = linspace(1,16,16);
h.TickLabels = {'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'};

title(h,'Soil Type')
zlim([0 800]);
print('-depsc', 'soiltype_1D_default')

%% --------1D SOIL TYPE DIFFERENCE --------

make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.05 0.045], [0.06 0.025], [0.1 0.2]);
if ~make_it_tight,  clear subplot;  end

figure(5)
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 450 600])

subplot(2,1,1)
surf(lon2d, lat2d, hgt,soilt_gssurgo); shading interp; view(2); 

ylabel('Latitude')
zlabel('Surface Elevation (m)')
ylim([35.1, 38.1])
xlim([-99.3,-95.6])
title_left('a) soil type from gSSURGO')

subplot(2,1,2)
surf(lon2d, lat2d, hgt,soilt_default); 
shading interp; view(2); 
xlabel('Longitude')
ylabel('Latitude')
zlabel('Surface Elevation (m)')
colormap (jet(16))
caxis([.5,16.5]) % value from the exact bottom to the exact top of the colorbar
zlim([0 800]);
ylim([35.1, 38.1])
xlim([-99.3,-95.6])
title_left('b) soil type from STATSGO')

hp3 = get(subplot(2,1,2),'Position');
h = colorbar('eastoutside')
set(h, 'Position', [hp3(3) + 0.125 hp3(2)+0.1 0.015  hp3(4)*1.5 ]);
% h.Label.String = ['SM pct change (%)'];
h.Ticks = linspace(1,16,16);
h.TickLabels = {'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'};

title(h,'Soil Type')

print('-dpdf', 'soiltype_1D_combined','-bestfit')
print('-depsc', 'soiltype_1D_combined')
%% ---------3D SOIL TYPE DATA ---------
figure(5)
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 600 600])
 
subplot(2,1,1)
surf(lon2d, lat2d, hgt,soilt_gssurgo); shading interp; view(20,50); 
 
% ylabel('Latitude')
zlabel('Surface Elevation (m)')
ylim([35.1, 38.1])
xlim([-99.3,-95.6])
title_left('a) soil type from gSSURGO')
 
subplot(2,1,2)
surf(lon2d, lat2d, hgt,soilt_default); 
shading interp; view(20,50);
xlabel('Longitude')
ylabel('Latitude')
zlabel('Surface Elevation (m)')
colormap (jet(16))
caxis([.5,16.5]) % value from the exact bottom to the exact top of the colorbar
zlim([0 800]);
ylim([35.1, 38.1])
xlim([-99.3,-95.6])
title_left('b) soil type from STATSGO')
 
hp3 = get(subplot(2,1,2),'Position');
h = colorbar('eastoutside')
set(h, 'Position', [hp3(3) + 0.15 hp3(2)+0.1 0.015  hp3(4)*1.5 ]);
% h.Label.String = ['SM pct change (%)'];
h.Ticks = linspace(1,16,16);
h.TickLabels = {'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'};
 
title(h,'Soil Type')
 
print('-dpdf', 'soiltype_3D_combined','-bestfit')
 
print('-depsc', 'soiltype_3D_combined')



%% Define 2d soil type matrix 16x16 with row = gssurgo, col = default 
labels = {'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'};

soilt_2d = nan(16,16);
for ii = 1:16
    for jj = 1:16 
        soilt_2d_count(ii,jj) = length(find (soilt_gssurgo == ii & soilt_default == jj));
    end 
end

soilt_pct(1,:) = sum(soilt_2d_count,1)/sum(soilt_2d_count(:))*100; % gssurgo
soilt_pct(2,:) = sum(soilt_2d_count,2)/sum(soilt_2d_count(:))*100; % statsgo
%%
figure
set(gcf,'PaperPositionMode', 'auto');
set(gcf,'Position', [100 100 600 250])

bar(1:16,soilt_pct');
set(gca, 'XTick', 1:16, ...                             % Change the axes tick marks
         'XTickLabel', labels, ...  %   and tick labels
         'TickLength', [0 0]);    
xtickangle(45)
grid on; 
legend('gSSURGO', 'STATSGO')
ylabel('percentage(%)')
print('-depsc', 'SOIL_TYPE_BARPLOT')

%%

labels = {'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'};
figure(10)
set(gcf,'PaperPositionMode', 'auto');
set(gcf,'Position', [100 100 600 600])
imagesc(soilt_2d_count); view(2);

colormap(flipud(gray)); 
grid on;
caxis([0 75000]);

textStrings = num2str(soilt_2d_count(:), '%d');       % Create strings from the matrix values
textStrings = strtrim(cellstr(textStrings));  % Remove any space padding
[x, y] = meshgrid(1:16);  % Create x and y coordinates for the strings
hStrings = text(x(:), y(:), textStrings(:), ...  % Plot the strings
                'HorizontalAlignment', 'center');
% set(hStrings, {'Color'}, num2cell(textColors, 2));  % Change the text colors

set(gca, 'XTick', 1:16, ...                             % Change the axes tick marks
         'XTickLabel', labels, ...  %   and tick labels
         'YTick', 1:16, ...
         'YTickLabel', labels, ...
         'TickLength', [0 0]);      
set(gca,'xaxisLocation','top')
xlabel('STAS');xtickangle(30)
ylabel('gSSURGO');ytickangle(20)
grid on;
set(gca, 'XTick', .5:15.5)
set(gca, 'YTick', .5:15.5)
print('-dpdf', 'gssurgo_default_soil_type')

%%
soilt_2d_pct = soilt_2d_count/sum(soilt_2d_count(:))*100;
figure(11)
set(gcf,'PaperPositionMode', 'auto');
set(gcf,'Position', [100 100 600 600])
imagesc(soilt_2d_pct); view(2);
grid on;
colormap(flipud(gray)); 

caxis([0 50]);

textStrings = num2str(soilt_2d_pct(:), '%0.2f');       % Create strings from the matrix values
textStrings = strtrim(cellstr(textStrings));  % Remove any space padding

for ii = 1:length(textStrings)
    if strcmp(textStrings(ii),'0.00')
        textStrings(ii) = {''};
    end 
end


[x, y] = meshgrid(1:16);  % Create x and y coordinates for the strings
hStrings = text(x(:), y(:), textStrings(:), ...  % Plot the strings
                'HorizontalAlignment', 'center');
% set(hStrings, {'Color'}, num2cell(textColors, 2));  % Change the text colors

set(gca, 'XTick', 1:16, ...                             % Change the axes tick marks
         'XTickLabel', labels, ...  %   and tick labels
         'YTick', 1:16, ...
         'YTickLabel', labels, ...
         'TickLength', [0 0]);      
set(gca,'xaxisLocation','top')
xlabel('STAS');xtickangle(30)
ylabel('gSSURGO');ytickangle(20)

grid on;
set(gca, 'XTick', .5:15.5)
set(gca, 'YTick', .5:15.5)
print('-dpdf', 'gssurgo_default_soil_type_pct')

%% calculate statistics
tot_default = sum(soilt_2d_pct);
tot_gssurgo = sum(soilt_2d_pct,2)';
for ii = 1:16
    diff(ii)    = tot_default(ii) - soilt_2d_pct(ii,ii);
    same(ii)    = soilt_2d_pct(ii,ii);
end
sum_same = sum(same);

%% 

diff_LH = ncread('./DATA_DIFF/diff_LH_1-9.nc', 'diff_LH');
diff_SH = ncread('./DATA_DIFF/diff_SH_1-9.nc', 'diff_SH');
diff_SM1 = ncread('./DATA_DIFF/diff_SM1_1-9.nc', 'diff_SOIL_M1');
diff_SM2 = ncread('./DATA_DIFF/diff_SM2_1-9.nc', 'diff_SOIL_M2');
diff_SM3 = ncread('./DATA_DIFF/diff_SM3_1-9.nc', 'diff_SOIL_M3');
diff_SM4 = ncread('./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_LH_soilt = nanmean(diff_LH(:,:,1,:), 4);
diff_SH_soilt = nanmean(diff_SH(:,:,1,:), 4);
diff_SM_soilt = nanmean(diff_SM(:,:,1,:), 4);

nanmean(diff_LH_soilt(:))
nanmean(diff_SH_soilt(:))
%  surf(lon2d, lat2d, diff_LH_soilt ); view(2); shading interp;

diff_LH_mean = zeros(16,16);
diff_SH_mean = zeros(16,16);

for ii = 1:16
    for jj = 1:16 
        diff_LH_mean(ii,jj) = nanmean(diff_LH_soilt (soilt_gssurgo == ii & soilt_default == jj));
        diff_SH_mean(ii,jj) = nanmean(diff_SH_soilt (soilt_gssurgo == ii & soilt_default == jj));
        diff_SM_mean(ii,jj) = nanmean(diff_SM_soilt (soilt_gssurgo == ii & soilt_default == jj));
    end 
end 
diff_LH_mean(isnan(diff_LH_mean)) = 0;
diff_SH_mean(isnan(diff_SH_mean)) = 0;
diff_SM_mean(isnan(diff_SM_mean)) = 0;

%%
labels = {'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'};

figure(12)
addpath('/Users/yang560/Dropbox/PNNL/Matlab_functions/utilities/')
set(gcf,'PaperPositionMode', 'auto');
set(gcf,'Position', [100 100 600 500])
imagesc(diff_LH_mean); view(2);
grid on;

h = colorbar
title(h,'LH diff. (W m^{-2})')
caxis([-60, 60]);colormap(bluewhitered); 

textStrings = num2str(diff_LH_mean(:), '%0.2f');       % Create strings from the matrix values
textStrings = strtrim(cellstr(textStrings));  % Remove any space padding

for ii = 1:length(textStrings)
    if strcmp(textStrings(ii),'0.00') || strcmp(textStrings(ii),'-0.00')
        textStrings(ii) = {''};
    end 
end


[x, y] = meshgrid(1:16);  % Create x and y coordinates for the strings
hStrings = text(x(:), y(:), textStrings(:), ...  % Plot the strings
                'HorizontalAlignment', 'center');
% set(hStrings, {'Color'}, num2cell(textColors, 2));  % Change the text colors

set(gca, 'XTick', 1:16, ...                             % Change the axes tick marks
         'XTickLabel', labels, ...  %   and tick labels
         'YTick', 1:16, ...
         'YTickLabel', labels, ...
         'TickLength', [0 0]);      
set(gca,'xaxisLocation','top')
xlabel('STAS');xtickangle(30)
ylabel('gSSURGO');ytickangle(20)

grid on;
set(gca, 'XTick', .5:15.5)
set(gca, 'YTick', .5:15.5)
print('-dpdf', 'LH_diff')


%%
weighted_diff_LH_mean = diff_LH_mean .* soilt_2d_pct/100;
weighted_diff_SH_mean = diff_SH_mean .* soilt_2d_pct/100;
weighted_diff_SM_mean = diff_SM_mean .* soilt_2d_pct/100;

make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.05 0.045], [0.06 0.025], [0.1 0.2]);
if ~make_it_tight,  clear subplot;  end

figure(13)
addpath('/Users/yang560/Dropbox/PNNL/Matlab_functions/utilities/')
set(gcf,'PaperPositionMode', 'auto');
set(gcf,'Position', [100 100 500 500])

imagesc(weighted_diff_LH_mean); view(2);

h=colorbar('southoutside')
title(h,'weighted LH diff. (W m^{-2})');
caxis([-1, 1]);colormap(bluewhitered); 

textStrings = num2str(weighted_diff_LH_mean(:), '%0.2f');       % Create strings from the matrix values
textStrings = strtrim(cellstr(textStrings));  % Remove any space padding

for ii = 1:length(textStrings)
    if strcmp(textStrings(ii),'0.00') || strcmp(textStrings(ii),'-0.00')
        textStrings(ii) = {''};
    end 
end


[x, y] = meshgrid(1:16);  % Create x and y coordinates for the strings
hStrings = text(x(:), y(:), textStrings(:), ...  % Plot the strings
                'HorizontalAlignment', 'center');
% set(hStrings, {'Color'}, num2cell(textColors, 2));  % Change the text colors

set(gca, 'XTick', 1:16, ...                             % Change the axes tick marks
         'XTickLabel', labels, ...  %   and tick labels
         'YTick', 1:16, ...
         'YTickLabel', labels, ...
         'TickLength', [0 0]);      
set(gca,'xaxisLocation','top')
xlabel('STATSGO');xtickangle(30)
ylabel('gSSURGO');ytickangle(20)

grid on;
set(gca, 'XTick', .5:15.5)
set(gca, 'YTick', .5:15.5)
print('-depsc', 'LH_diff_weighted')

%%
figure (14)
addpath('/Users/yang560/Dropbox/PNNL/Matlab_functions/utilities/')
set(gcf,'PaperPositionMode', 'auto');
set(gcf,'Position', [100 100 500 500])
imagesc(weighted_diff_SH_mean); view(2);
h = colorbar('southoutside')
title(h,'weighted SH diff. (W m^{-2})');
caxis([-1, 1]);colormap(bluewhitered); 

textStrings = num2str(weighted_diff_SH_mean(:), '%0.2f');       % Create strings from the matrix values
textStrings = strtrim(cellstr(textStrings));  % Remove any space padding

for ii = 1:length(textStrings)
    if strcmp(textStrings(ii),'0.00') || strcmp(textStrings(ii),'-0.00')
        textStrings(ii) = {''};
    end 
end


[x, y] = meshgrid(1:16);  % Create x and y coordinates for the strings
hStrings = text(x(:), y(:), textStrings(:), ...  % Plot the strings
                'HorizontalAlignment', 'center');
% set(hStrings, {'Color'}, num2cell(textColors, 2));  % Change the text colors

set(gca, 'XTick', 1:16, ...                             % Change the axes tick marks
         'XTickLabel', labels, ...  %   and tick labels
         'YTick', 1:16, ...
         'YTickLabel', labels, ...
         'TickLength', [0 0]);      
set(gca,'xaxisLocation','top')
xlabel('STATSGO');xtickangle(30)
ylabel('gSSURGO');ytickangle(20)

grid on;
set(gca, 'XTick', .5:15.5)
set(gca, 'YTick', .5:15.5)
 print('-depsc', 'SH_diff_weighted')

%% 
figure (15)
addpath('/Users/yang560/Dropbox/PNNL/Matlab_functions/utilities/')
set(gcf,'PaperPositionMode', 'auto');
set(gcf,'Position', [100 100 500 500])
imagesc(diff_SM_mean); view(2);
h=colorbar('southoutside')
title(h,' SM diff. (m)');
caxis([-.8, .8]);colormap(bluewhitered); 

textStrings = num2str(diff_SM_mean(:), '%0.2f');       % Create strings from the matrix values
textStrings = strtrim(cellstr(textStrings));  % Remove any space padding

for ii = 1:length(textStrings)
    if strcmp(textStrings(ii),'0.00') || strcmp(textStrings(ii),'-0.00')
        textStrings(ii) = {''};
    end 
end


[x, y] = meshgrid(1:16);  % Create x and y coordinates for the strings
hStrings = text(x(:), y(:), textStrings(:), ...  % Plot the strings
                'HorizontalAlignment', 'center');
% set(hStrings, {'Color'}, num2cell(textColors, 2));  % Change the text colors

set(gca, 'XTick', 1.:16, ...                             % Change the axes tick marks
         'XTickLabel', labels, ...  %   and tick labels
         'YTick', 1:16, ...
         'YTickLabel', labels, ...
         'TickLength', [0 0]);      
set(gca,'xaxisLocation','top')
xlabel('STATSGO');xtickangle(30)
ylabel('gSSURGO');ytickangle(20)
% 
% for ii = 1:16    
%         plot(1:16, ii*ones(1,16),':');hold on; 
%         plot(ii*ones(1,16), 1:16, ':');hold on;
% end


grid on; 
set(gca, 'XTick', .5:15.5)
set(gca, 'YTick', .5:15.5)

print('-depsc', 'SM_diff')

 
 
%%
sum(weighted_diff_LH_mean(:))
sum(weighted_diff_SH_mean(:))

for ii = 1:16
    for jj = 1:16 
        if ii == jj 
            weight_diff_LH_mean(ii,jj) = 0;
        end 
    end 
end
sum(weighted_diff_LH_mean(:))

sum(weighted_diff_SH_mean(:))


%% 
figure(16)
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 600 300])
surf(lon2d, lat2d, hgt,diff_LH_soilt); 
shading interp; view(20,50);
xlabel('Longitude')
ylabel('Latitude')
zlabel('Surface Elevation (m)')
caxis([-30,30]);colormap (bluewhitered);
h = colorbar;
 % value from the exact bottom to the exact top of the colorbar

% h.Ticks = linspace(1,16,16);
% h.TickLabels = {'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'};

title(h,'LH diff. (W m^{-2})')
zlim([0 800]);
print('-depsc', 'diff_LH_3D_by_gssurgo')

nanmean(diff_LH_soilt(:))

%%
figure(17)
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 600 400])
surf(lon2d, lat2d, hgt,diff_LH_soilt); 
shading interp; view(20,50);view(2);
xlabel('Longitude')
ylabel('Latitude')
zlabel('Surface Elevation (m)')
caxis([-30,30]);colormap (bluewhitered);
h = colorbar;
 % value from the exact bottom to the exact top of the colorbar

% h.Ticks = linspace(1,16,16);
% h.TickLabels = {'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'};

title(h,'LH diff. (W m^{-2})')
zlim([0 800]);
print('-depsc', 'diff_LH_2D_by_gssurgo')


%% 
diff_LH_soilt_sand = diff_LH_soilt;
diff_LH_soilt_sand(~(soilt_gssurgo ~= 1 & soilt_default == 1)) = 0;
figure(18)
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 600 400])
surf(lon2d, lat2d, hgt,diff_LH_soilt_sand); 
shading interp; view(20,50);view(2);
xlabel('Longitude')
ylabel('Latitude')
zlabel('Surface Elevation (m)')
caxis([-30,30]);colormap (bluewhitered);
h = colorbar;

title(h,'LH diff. (W m^{-2})')
print('-dpdf', 'diff_LH_2D_sand')

%%
diff_LH_soilt_sand = diff_LH_soilt;
diff_LH_soilt_sand(~(soilt_gssurgo == 15 & soilt_default ~= 15)) = 0;
figure(19)
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 600 400])
surf(lon2d, lat2d, hgt,diff_LH_soilt_sand); 
shading interp; view(20,50);view(2);
xlabel('Longitude')
ylabel('Latitude')
zlabel('Surface Elevation (m)')
caxis([-30,30]);colormap (bluewhitered);
h = colorbar;

title(h,'LH diff. (W m^{-2})')
print('-dpdf', 'diff_LH_2D_bedrock')


