
% Script to apply an empirically derived depth-reflectance scaling
% relationship to derive depths for Sentinel-2. Here, the scaling
% relationship was derived by comparing the Landsat-8-derived lake depths
% (using the physically-based methods of Pope et al., 2016) with the
% Sentinel-2 red band surface reflectance values and then applying a
% regression relationship to these data. Various regressions were tested
% with the surface reflectance data (for green, red and blue bands), but we
% found that a 1-term power law regression applied to the red band data
% produced the best statistical relationship. 

% Power law regression is y = a * x^-b 
% Here, a = 0.2764 and b = -0.8952 
% Regression values are y = (1.0528 * 10^3) * (x)^-0.8952 
% where depth is y and x is the red band surface reflectance (in digital
% numbers, not scaled by scalar of 10,000) 

% Note that this script takes a LONG time to run. 

% Andrew Williamson, University of Cambridge, September 2017 

%% Specify input directories and parameters 

main_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/'; 
S2_folder = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_S2';
S2_bands_folder = [S2_folder '/bands'];
S2_NDWI_folder = [S2_folder '/NDWI_masks'];
S2_mask_folder = [main_dir '/masks']; 
S2_out_depth_folder = [S2_folder '/depths']; 
S2_out_info_folder = [S2_folder '/lake_info']; 
S2_pix_siz = 0.01; S2_pix_siz_m = 10; % specifying S2 resolution  

%% Make directories for saving files 

mkdir(S2_out_depth_folder); mkdir(S2_out_info_folder);  

%% Load data into workspace 

cd(S2_bands_folder); 
S2_red_list = dir('*_red.mat'); 
cd(S2_NDWI_folder); 
S2_NDWI_list = dir('*_NDWI_SGL_mask.mat'); 
cd(S2_mask_folder); 
GrIS_mask = imread('S2_mask_10m.tif'); 

num_img = length(S2_red_list); 
string = num2str(num_img); 
disp(['Calculating lake depths for ' string ' Sentinel-2 images.'])
disp('Lake depths being calculated with empirically derived depth-reflectance relationship.') 

for i = 1:num_img 
    
    % get data and write image details to command line 
    dates = {S2_red_list(:).name}; 
    dates = dates'; % transpose matrix as it's nicer this way
    current_date = dates{i,1}; 
    year = current_date(1,4:7);
    month = current_date(1,8:9); 
    day = current_date(1,10:11); 
    curr_date_str = [day '/' month '/' year]; % format date string nicely 
    curr_imnum_str = num2str(i); 
    disp(' ')
    disp(['Now analysing image from ' curr_date_str ', image number '...
        curr_imnum_str ' of ' string '.'])
    
    % new variable for short name string for use later 
    short_dt_str = [year month day]; 
    
    % get the filenames for loading reflectance data and masks 
    red_fname = S2_red_list(i).name; 
    NDWI_fname = S2_NDWI_list(i).name; 
    
    % load the data into the workspace (note that the pan_band and red_band
    % and NDWI_defined_SGLs names below are based on how the original
    % variables were saved to the .mat files - change if needed) 
    cd(S2_bands_folder); 
    disp('Loading red band surface reflectance data.')
    load(red_fname); 
    red_refl = red_band; 
    red_refl = red_refl ./ 10000; % scalar for S2 is 10,000 
    clear red_band 
    cd(S2_NDWI_folder); 
    disp('Loading NDWI-derived lake area mask.')
    load(NDWI_fname);
    lake_mask = NDWI_defined_SGLs;  
    
    % convert surface reflectance to double format for analysis and mask
    % no-data values (zeroes in S2), also mask the GrIS edge  
    red_refl = double(red_refl); 
    red_refl(red_refl == 0) = NaN; 
    red_refl(GrIS_mask == 0) = NaN; 
    lake_mask(GrIS_mask == 0) = NaN; 
    
    % get a lake-number-labelled NDWI mask 
    lake_mask(isnan(lake_mask)) = 0; 
    [lake_mask_labelled,lakenum] = bwlabel(lake_mask,8); % default 8 connectivity 
    
    % do the lake depth calculations for each lake in the array
    disp(['Calculating depths for all ' num2str(lakenum) ' lakes on the image.']);
    
    red_refl2 = red_refl; 
    red_refl2(lake_mask == 0) = NaN; 
    a = 0.2764; 
    tmpdepth_red2 = red_refl2 .^-0.8952; 
    tmpdepth_red1 = lake_mask .* a; 
    tmp_depth = tmpdepth_red2 .* tmpdepth_red1; 
    lake_depths = tmp_depth; 
    
    % mask the lake depth arrays off the ice-sheet edge with NaNs 
    lake_depths(GrIS_mask == 0) = NaN; 
    
    % get new variable names from above 
    new_depth_var_name = ['S2_' short_dt_str '_depths'];
    %new_info_var_name = ['S2_' short_dt_str '_lake_info'];
    
    % save the output files 
    cd(S2_out_depth_folder); 
    save([new_depth_var_name '.mat'],'lake_depths')
%     cd(S2_out_info_folder); 
%     save([new_info_var_name '.mat'],'lakeinfo')
%     
    % display progress update to command line 
    disp(' '); 
    progress = 100/num_img * i; % percentage based on number of images 
    
    figure(1); 
    imagesc(lake_depths) 
    
    disp(['Lake depth calculations now ' num2str(progress,'%.1f') '% complete.']) 
    
end 

cd(main_dir) 
