
% Script to compare Landsat-8 depths derived using a physically based
% method applied to the red and panchromatic band surface reflectance data
% against Sentinel-2 depths derived using the same physically based method
% applied to the Sentinel-2 red band. 

% Note that the script only does these comparisons for overlapping imagery
% - it is NOT designed to calculate Sentinel-2 lake depths for every image
% in a season; it is purely for validation of the Sentinel-2 method.

% Script has the option either to calculate Sentinel-2 depths within it or,
% if these calculations have already been done previously, to simply load
% the existing arrays, to save (a lot of) time. 

% Script assumes the Landsat-8 depths have already been calculated for
% these two images using the L8_lake_depth.m script, which is based on the
% physically based methods of Pope et al. (2016). These depths are to be
% regarded as ground-truth data (Williamson et al., 2017). Ultimately, once
% this validation script has been used, the calculations for Sentinel-2
% depth may be used in an entirely new script (S2_lake_depth_physical.m) to
% calculate Sentinel-2 lake depths for every image in a season. 

% Andrew Williamson, University of Cambridge, September 2017 

%% Specify image number, etc. 

pre_process = 'Y'; % do you need to pre-process the data? 'Y' if the S2 
% depth data have not already been calculated or 'N' if thay have been and
% can just be loaded from the directory; note that doing the calculations
% takes a lot of time, so always ensure this option is checked only if
% required! 
num_dates = 2; % can set up to do for 1 or 2 dates simultaneously 
date1 = '20160731'; % for the lake depth and surface reflectance band data
date2 = '20160701'; % for the lake depth and surface reflectance band data
depth_interp_type = 'nearest'; % 'nearest' or 'bilinear' for testing 
radius = 200; % amount erode ice edge by; a scalar to be multiplied by 10 m
S2_scalar = 10000; % for converting digital numbers to TOA reflectance 
root_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/'; 
testing_dir = [root_dir 'outputs_S2_L8_dpth_comp/']; mkdir(testing_dir); 
S2_pix_siz = 0.01; S2_pix_siz_m = 10; % specifying S2 resolution  
g_red = 0.8304; % g value for S2 red band calculated with methods of Pope
% and others (2016)
close all % close all open figure windows 
clc % clear command line 

%% Load data if not already done 

% specify at the command line the number of dates being analysed 
if num_dates == 1 
    loop = 1;
    disp('There is one date to analyse.') 
elseif num_dates == 2 
    loop = 2; 
    disp('There are two dates to analyse.') 
end
  
% depending on whether doing 1 or 2 dates, this will change to 1 or 2
% loops... 
for loop_occs = 1:loop 
    
if loop_occs == 1
    date = date1; 
elseif loop_occs == 2 
    date = date2; 
end
    
disp(['Doing analysis for image from ' date '.'])  

% Landsat-8 lake depth data 
load([root_dir 'outputs_LS8/depths/L8_' date '_depths.mat']);
S2_size = [10980 10980]; % Sentinel-2 image size 
lake_depths_10m = imresize(lake_depths,S2_size,depth_interp_type);
disp('Landsat-8 lake depth data loaded and interpolated for this date.')

% lake area masks from NDWI analysis 
load([root_dir 'outputs_S2/NDWI_masks/S2_' date '_NDWI_SGL_mask.mat']);
load([root_dir 'outputs_LS8/NDWI_10m_masks/L8_' date '_NDWI_10m_SGL_mask.mat']);
L8_lmask = L8_lmask_10m; 
S2_lmask = NDWI_defined_SGLs; 
disp('NDWI lake masks from Sentinel-2 and Landsat-8 loaded for this date.') 

%% Masking 

% do masking according to ice-sheet mask so that don't get 0 depth values
% from Landsat 8 
% must re-load the mask every time to avoid it gradually being eroded each
% time the script is run... 
disp('Masking ice-marginal areas using GIMP mask (Howat et al., 2014).') 
clear S2_mask_10m 
S2_mask_10m = imread([root_dir 'masks/S2_mask_10m.tif']);
e = strel('disk',radius); % to get rid of the issues at the ice edge 
S2_mask_10m = imerode(S2_mask_10m,e); % erode mask by radius value 

% mask these values as no-data values 
L8_lmask(S2_mask_10m == 0) = NaN; 
S2_lmask(S2_mask_10m == 0) = NaN;

% specify erosion amount at command line 
erosion_amt = (radius * 10) / 1000; % multiply by 10 m then get to km 
disp(['Ice sheet masking done, including erosion of ice edge by ' ...
    num2str(erosion_amt) ' km.'])

%% Do the Sentinel-2 lake depth calculations with the physical method 
if strcmp(pre_process,'Y') == 1

disp('Sentinel-2 data NOT pre-processed, so doing now.') 
% load the surface reflectance data 
load([root_dir 'outputs_S2/bands/S2_' date '_red.mat']);
disp('Sentinel-2 red band TOA data loaded for lake depth calculations.')
red_band = double(red_band); % convert to double format if not already 
red_band = red_band ./ S2_scalar; % scale data accordingly

% get darkest pixel value from the red band 
rinf_red = min(min(red_band)); 
rinf_red_str = num2str(rinf_red,'%.3f'); 
disp(['Darkest red band pixel value is ' rinf_red_str '.']); 

if rinf_red > 0.1 
    disp('Warning: darkest red band value is high - manually check.'); 
end 

% set up lake mask array 
lake_mask = S2_lmask; 

% get a lake-number-labelled NDWI mask 
lake_mask(isnan(lake_mask)) = 0; 
[lake_mask_labelled,lakenum] = bwlabel(lake_mask,8); % default 8 connectivity 
    
% set up blank arrays for depth information 
target_size = size(red_band); 
lake_depths = zeros(target_size); 
lakeinfo = zeros(lakenum,8); 
e = strel('disk',3); % to dilate lake mask by 3 pixels (bigger than the
% usual 1, to ensure that we are not including the shallowest water that
% has been omitted from the lake perimeters) 

% do the lake depth calculations for each lake in the array
disp(['Calculating depths for all of the ' num2str(lakenum) ' lakes on the image.']);
for n = 1:lakenum % defines number of times to loop 
        
    % show the lake number being worked on at the command line 
    disp(['Now working on lake number ' num2str(n) ' of a total ' num2str(lakenum) ' lakes.'])

    lmask = zeros(target_size); % creates lake mask with size as above  
    tmpref_red = red_band; % tmpref is equal to values of red band 
    
    % Get the individual lake mask
    lmask(lake_mask_labelled == n) = 1;
    
    % put lake area information into the lakeinfo array 
    lakeinfo(n,1) = n; % # of lake to 1st col of 'lakeinfo' array
    lakeinfo(n,2) = sum(sum(lmask)); % sum of pixels within lake
    lakeinfo(n,7) = (sum(sum(lmask)) * S2_pix_siz^2); % lake area in km^2

    % Dilate lake mask by one ring of pixels to find bare-ice albedo 
    lmaskdil = imdilate(lmask,e); % e is defined above - structuring element
    ledge = lmaskdil - lmask; % just get the edge pixels 
    edgepix_red = (red_band(ledge == 1)); % find the red edge pixel refl 
    ad_red = nanmean(edgepix_red); % take the mean of these values (ignore NaNs) 
    lakeinfo(n,5) = ad_red; % put this in lake info file 
    
    % mask reflectance arrays 
    tmpref_red(lmask == 0) = 0;  
   
    % calculate lake depths with methods of Sneed and Hamilton (2007); 
    % see details in Pope and others (2016) or Williamson and others (2017)  
    tmpdepth_red = (log(ad_red - rinf_red) - log(tmpref_red - rinf_red)) / (-g_red);   
    
    % add the two depth values together and take the average of them 
    tmpdepth = tmpdepth_red;   
    tmpdepth(lmask == 0) = 0; % if lmask value = 0, then tmpdepth = 0
    tmpdepth = -tmpdepth; % negative of tmpdepth to correct for -g
    lakeinfo(n,3) = sum(sum(tmpdepth)); % sum of all pixel depths
    lakeinfo(n,4) = max(max(tmpdepth)); % max lake depth 
    lakeinfo(n,6) = mean(mean(tmpdepth)); % mean lake depth 
    lakeinfo(n,8) = (sum(sum(tmpdepth)) * (S2_pix_siz_m)^2); % volume (m^3)
    
    % add calculated depths to lake depth array
    lake_depths = lake_depths + tmpdepth;
    
end % end of lake-depth calculation loop  
    
% mask the lake depth arrays off the ice-sheet edge with NaNs 
lake_depths(S2_mask_10m == 0) = NaN; 
    
% get new variable names from above 
new_depth_var_name = ['S2_' date '_depths'];
new_info_var_name = ['S2_' date '_lake_info'];

% save the output files 
cd(testing_dir)
save([new_depth_var_name '.mat'],'lake_depths')
cd(testing_dir); 
save([new_info_var_name '.mat'],'lakeinfo')

elseif strcmp(pre_process,'N') == 1 

disp('Sentinel-2 data have been pre-processed, so simply loading the existing array.') 
cd(testing_dir); 
load([testing_dir 'S2_' date '_depths.mat']);
S2_depth = lake_depths; 
    
end % end of pre-processing loop 

%% Identify locations for analysis 

% get locations that are identified as water cover both in S2 and L8 
% since the individual masks are binary arrays, adding together produces 2s
% where they overlap 
summed_lmask = L8_lmask + S2_lmask; 

% set up a comparisons array with size of 1 row and 2 columns 
S2_L8_comps = zeros(1,2); 
comps_row = 0; % set this up as a 0 to be changed within the for-loop 
L8_col = 1; % column to put the L8 lake depth in 
S2_col = 2; % column to put the S2 lake depth in 
%row_identifier = 3; % row identification check if wanted 
%col_identifier = 4; % column identification check if wanted 

%% Do the main steps of the analysis 

% for each row of the summed lake mask... 
for i = 1:size(summed_lmask,1)
    for j = 1:size(summed_lmask,2) % and for each column... 
        
        % check the summed_lmask to find the overlap identification (of 2) 
        overlap_identifier = summed_lmask(i,j); 
        
        % if this is 2 (i.e. both S2 and L8 data exist) 
        if overlap_identifier == 2 
            
            % get the L8 lake depth (in metres) from the right row and
            % column (i,j) 
            L8_depth_val = lake_depths_10m(i,j); 
            
            % get the S2 lake depth (in metres) from the right row and
            % column (i,j) 
            S2_depth = lake_depths; 
            S2_depth_val = S2_depth(i,j); 

            % add 1 to make a new row in this array 
            comps_row = comps_row + 1; 
            % put L8 lake depth in metres the 1st col  
            S2_L8_comps(comps_row,L8_col) = L8_depth_val; 
            % put S2 reflectance value (raw digital numbers) into 2nd col
            S2_L8_comps(comps_row,S2_col) = S2_depth_val; 
            
            % included in case manual checks wanted but now not needed 
            % S2_L8_comps(comps_row,row_identifier) = i; 
            % S2_L8_comps(comps_row,col_identifier) = j; 
        
        end % end of if loop 
    end % end of j loop 
end % end of i loop 

% test to state that if this is the first loop ('loop_occs') or the second
% loop, and the number of loops to do is 2, then name the array
% appropriately 
% if the loop number is only 1 and there is only 1 loop to do, then set it
% just as 1 
if loop_occs == 2 && loop == 2 
    S2_L8_comps_date2 = S2_L8_comps;  
elseif loop_occs == 1 && loop == 2 
    S2_L8_comps_date1 = S2_L8_comps;  
elseif loop_occs == 1 && loop == 1 
    S2_L8_comps_date1 = S2_L8_comps; 
end

end 

%% Make scatter plot 

disp('Making scatter plot of data.')

ft_size = 18; % specify font size for entire figure 

hFig = figure('rend','painters','pos',[10 10 1400 1400]); % new figure  
% original plotting without transparency (data not visible) 
%data = scatter(x_variable,y_variable,'MarkerEdgeColor','r','Marker','.'); 
set(hFig,'Name','Sentinel-Landsat lake depth comparison','PaperType',...
    'a3','PaperOrientation','landscape','PaperUnits','normalized')
transparency_level = 0.05;

if num_dates == 1 
    num_plot_times = 1; 
elseif num_dates == 2 
    num_plot_times = 2; 
end

for i = 1:num_plot_times 
  
    if i == 1 
    % copy this data to leave the original unchanged
        scatter_data = S2_L8_comps_date1; 
    elseif i == 2 
        hold on 
        scatter_data = S2_L8_comps_date2;
    end

% clean the data to remove where L8 predicted negative values or 0s 
scatter_data(scatter_data(:,1) <= 0,:) = [];
scatter_data(scatter_data(:,2) <= 0, :) = []; 
scatter_data(isnan(scatter_data),:) = []; 
scatter_data(isinf(scatter_data),:) = []; 

% get tidy scatter data 
y_variable = scatter_data(:,L8_col); % extract the x variable 
x_variable = scatter_data(:,S2_col); % extract the y variable 
n = size(scatter_data,1); % find number of data points 

if num_plot_times == 2 && i == 1 
    x_variable1 = x_variable; 
    y_variable1 = y_variable; 
elseif num_plot_times == 2 && i == 2
    x_variable2 = x_variable; 
    y_variable2 = y_variable; 
end 

if i == 1 
    colr = 'r'; 
    type = 's';
elseif i == 2 
    colr = 'b'; 
    type = 'o'; 
end 

data = scatter(x_variable,y_variable,15,type,'MarkerFaceAlpha',...
    transparency_level,'MarkerFaceColor',colr,'MarkerEdgeColor','none');
ylabel('Landsat 8 lake depth (m)','FontSize',ft_size); 
y = xlabel('Sentinel-2 lake depth (m)','FontSize',ft_size); 
%set(y, 'position', get(y,'position')-[0.1,0,0]); % offset y label 

end 

%% Sort out the regression data first 

if num_plot_times == 1
    reg_x_variable = x_variable; 
    reg_y_variable = y_variable; 
elseif num_plot_times == 2
    reg_x_variable = [x_variable1;x_variable2]; 
    reg_y_variable = [y_variable1;y_variable2];
end 

%% Add power regression relationship to the graph 

hold on 

reg_to_add = 'polynomial'; % which regression relationship do you want? either
% 'power, 'polynomial' (which is 1st order) or 'both' 

% power law curve 
[f,gof,output] = fit(reg_x_variable,reg_y_variable,'power1');
a = f.a; % a value from power law relation for the data y = a * x^b; 
b = f.b; % b value from power law relation for the data y = a * x^b; 

% add power law curve to the data 
xl = xlim; 
x_low = xl(1,1); % get the x-axis lower limit
x_high = xl(1,2); % get the x-axis higher limit 
x_locs = linspace(x_low,x_high,1000); % get the plotting locations for x
y_pred = a * (x_locs).^(b); % predict y values for these x locations 

if strcmp(reg_to_add,'power') == 1 || strcmp(reg_to_add,'both') == 1
    plot(x_locs,y_pred,'Color','k','LineWidth',1); 
end

ylim([0 6]) % set to the limits of the depth data 
xlim([0 6]) % set x limits the same 

% extract the stats from goodness of fit data 
r2 = gof.rsquare; 
df = gof.dfe; 
rmse = gof.rmse; 

% add text in a legend to show the p-value and r-value 
str = {['R^2 = ' num2str(r2,'%.3f')],['RMSE = ' num2str(rmse,'%.3f') ' m'],...
    ['df = ' num2str(df)]};
dim = [0.765 0.6 0.28 0.3]; % specify the dimensions and location of text box 
annotation('textbox',dim,'String',str,'FitBoxToText','on',...
    'BackgroundColor','w','FontSize',ft_size); 

% change to 1 dp on the y-axis for Landsat-8 depths 
ax = gca; 
ax.YAxis.TickLabelFormat = '%,.2f'; 
ax.XAxis.TickLabelFormat = '%,.2f';

% make font size of font size and put x- and y-grids on 
set(gca,'XGrid','on','YGrid','on','FontSize',ft_size)
set(gca,'TickLength',[0.005 0.005]); % tick length 
box on % box around figure 

%% Add linear regression = 1st order polynomial 

% Pearson linear regression 
[r,p] = corr(reg_x_variable,reg_y_variable,'Type','Pearson','rows','pairwise');

% get the polynomial curve values for the regression variables 
polyf = polyfit(reg_x_variable,reg_y_variable,1); 

% add polynomial curve to graph 
xl = xlim; 
x_low = xl(1,1); % get the x-axis lower limit
x_high = xl(1,2); % get the x-axis higher limit 
x_locs = linspace(x_low,x_high,1000); % get the plotting locations for x
y1 = polyval(polyf,x_locs); 

if strcmp(reg_to_add,'polynomial') == 1 || strcmp(reg_to_add,'both') == 1
    plot(x_locs,y1,'Color','k','LineWidth',1); 
end

%% Save figure 

disp('Saving figure.') 
filename = 'S2_L8_physical_depth_comp'; 
cd(testing_dir);
saveas(gcf,[filename '.png']); 
cd(root_dir); 
