clc 
clear 
close all 

%% PLOT the 3D soil moisture on surface. 
addpath('../ncfiles_10km/')
addpath('/Users/yang560/Dropbox/PNNL/Matlab_functions/utilities/')

soilt_gssurgo = ncread('geo_em.d01.nc.10km', 'SCT_DOM');
lat2d     = ncread('geo_em.d01.nc.10km', 'XLAT_M');
lon2d     = ncread('geo_em.d01.nc.10km','XLONG_M');
hgt       = ncread('geo_em.d01.nc.10km', 'HGT_M');

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

%---- add  soil moisture file ----
soil_1 = ncread('seavar_SM1_19-23.nc','seavar_SOIL_M1');
soil_2 = ncread('seavar_SM2_19-23.nc','seavar_SOIL_M2');
soil_3 = ncread('seavar_SM3_19-23.nc','seavar_SOIL_M3');
soil_4 = ncread('seavar_SM4_19-23.nc','seavar_SOIL_M4');

soil_1(:,:,:,5) = squeeze(mean(soil_1 , 4));
soil_2(:,:,:,5) = squeeze(mean(soil_2 , 4));
soil_3(:,:,:,5) = squeeze(mean(soil_3 , 4));
soil_4(:,:,:,5) = squeeze(mean(soil_4 , 4));

soil_1_diff = func_calc_diff(soil_1);
soil_2_diff = func_calc_diff(soil_2);
soil_3_diff = func_calc_diff(soil_3);
soil_4_diff = func_calc_diff(soil_4);

%% ==== saveas netcdf ====
func_saveasnetcdf(lat2d,lon2d,soil_1_diff,'soil_1_diff');
func_saveasnetcdf(lat2d,lon2d,soil_2_diff,'soil_2_diff');
func_saveasnetcdf(lat2d,lon2d,soil_3_diff,'soil_3_diff');
func_saveasnetcdf(lat2d,lon2d,soil_4_diff,'soil_4_diff');


%% ----- PLOT SECTION -----
% close all
% func_plot_soil(lon2d, lat2d, hgt,soil_1,'SM1'); print('-depsc', 'SM1_4s') 
% func_plot_soil(lon2d, lat2d, hgt,soil_2,'SM2'); print('-depsc', 'SM2_4s')
% func_plot_soil(lon2d, lat2d, hgt,soil_3,'SM3'); print('-depsc', 'SM3_4s')
% func_plot_soil(lon2d, lat2d, hgt,soil_4,'SM4'); print('-depsc', 'SM4_4s')

%% related to resolution impact, absolute difference
close all
season = 5;
func_plot_soil_res(lon2d, lat2d, hgt,soil_1_diff,...
                                     soil_2_diff,...
                                     soil_3_diff,...
                                     soil_4_diff,season)

%% related to resolution impact, but in terms of percentage

[sm1_mean_pct_change] = func_cal_domain_mean(soil_1, soil_1_diff);
[sm2_mean_pct_change] = func_cal_domain_mean(soil_2, soil_2_diff);
[sm3_mean_pct_change] = func_cal_domain_mean(soil_3, soil_3_diff);
[sm4_mean_pct_change] = func_cal_domain_mean(soil_4, soil_4_diff);

func_plot_soil_pct_change_res(lon2d, lat2d, hgt,sm1_mean_pct_change,sm2_mean_pct_change,sm3_mean_pct_change,sm4_mean_pct_change);

% %%
% % close all
% % season = 5;
% % func_plot_soil_vert_profile(lon2d, lat2d, hgt,soil_1_diff,...
% %                                                   soil_2_diff,...
% %                                                   soil_3_diff,...
% %                                                   soil_4_diff,season)
% %%
% 
% close all;
% season = 5;
% func_plot_soil_pdf(soil_1_diff,soil_1,...
%                    soil_2_diff,soil_2,...
%                    soil_3_diff,soil_3,... 
%                    soil_4_diff,soil_4,season)
% %%
% close all 
[aspect,slope,gradN,gradE] = gradientm(lat2d,lon2d,hgt);
% 
% func_plot_soil_scatter(soil_1_diff,soil_2_diff,soil_3_diff,soil_4_diff,slope);
% print('-depsc', '10km_topo_slope')
%                               
% %% Calculate and plot the bar plot for each slope bin
% func_plot_soil_slope_bar(soil_1_diff,soil_2_diff,soil_3_diff,soil_4_diff,slope);
% 
% print('-depsc', '10km_slope_barplot')


%%
func_plot_soil_slope_bar_p3(soil_1_diff,soil_2_diff,soil_3_diff,soil_4_diff,slope);
print('-depsc', '10km_slope_barplot_p3')