
% Script to compare Landsat-8 depths against surface reflectance data from
% Sentinel-2 to generate a scatter plot to show the data with a power law
% regression added. 

% Possible to run for a variety of different bands (red, blue, green) and 
% for different imagery dates assuming there are image pairs (only works
% with two atm). 

% N.B. script only deals with power law or linear regression; anything more
% complex that may be established is done with the data eleswhere 
% (e.g. SPSS). 

% Andrew Williamson, University of Cambridge, September 2017 

%% Specify image number etc. 

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
data_loaded = 'N'; % to save time if data already in workspace - need to 
% leave this as 'N' all the time if doing 2 dates! 
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; 
band_tp = 'red'; % band to do the analysis for; either 'blue','red' or 'green'
root_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/'; 
out_dir = [root_dir 'output_figs/']; 
mkdir(out_dir); 
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 ' for ' band_tp ' band.']) 

if strcmp(data_loaded,'N') == 1 

    % lake depth data 
    load([root_dir 'outputs_LS8/depths/L8_' date '_depths.mat']);
    S2_size = [10980 10980]; 
    lake_depths_10m = imresize(lake_depths,S2_size,depth_interp_type);
    disp('Landsat-8 lake depth data loaded and interpolated.')
    
    % surface reflectance data 
    load([root_dir 'outputs_S2/bands/S2_' date '_blue.mat']);
    load([root_dir 'outputs_S2/bands/S2_' date '_green.mat']);
    load([root_dir 'outputs_S2/bands/S2_' date '_red.mat']);
    disp('Sentinel-2 surface reflectance data loaded.')
    
    % 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.') 
    
elseif strcmp(data_loaded,'Y') == 1
    
    disp('Depth and band data already in workspace, so moving on to analysis.') 

end 

red_blue_band = red_band + blue_band; 
red_band = red_band ./ S2_scalar;
green_band = green_band ./ S2_scalar;
blue_band = blue_band ./ S2_scalar;
red_blue_band = red_blue_band ./ S2_scalar;

%% 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... 
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.'])

%% 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); 
            
            % depending on which band you want to do, then get the S2
            % reflectance value from the appropriate band for the same
            % row and column (i,j) 
            if strcmp(band_tp,'red') == 1 
                
                S2_refl_val = red_band(i,j); 
                
            elseif strcmp(band_tp,'green') == 1 
                 
                S2_refl_val = green_band(i,j); 
                
            elseif strcmp(band_tp,'blue') == 1 
                
                S2_refl_val = blue_band(i,j); 
                
            elseif strcmp(band_tp,'red_blue') == 1 
                
                S2_refl_val = red_blue_band(i,j); 
                
            end
            
            % 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_refl_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 regression','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(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 ' band_tp ' band TOA reflectance'],...
    'FontSize',ft_size); 
%set(y, 'position', get(y,'position')-[0.1,0,0]); % offset y label 

end 

%% Add regression relationship to the graph 

hold on 

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 
    
% Pearson linear regression 
% [r,p] = corr(reg_x_variable,reg_y_variable,'Type','Pearson','rows','pairwise');

% 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 
plot(x_locs,y_pred,'Color','k','LineWidth',1); 

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

% 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 = '%,.1f'; 
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 

%% Save figure 

disp('Saving figure.') 
filename = ['S2_L8_regression_' band_tp]; 
cd(out_dir);
saveas(gcf,[filename '.png']); 
cd(root_dir); 
