% Script to calculate lake depths and volumes for a series of Landsat-8
% images and to save them. 

% Depths calculated using the methods of Pope and others (2016), i.e. by
% taking the average of the red- and panchromatic-band-derived depths to
% calculate an individual pixel depth. Original function based on that in
% use by Banwell and others (2014) and Arnold and others (2014), also in
% use by Williamson and others (2017). 

% Points to the folders with the bands and NDWI masks within the script.
% Within these folders, the files need to be saved as follows: 
% L8_20160604_red.mat for red band (exact same date format, etc.) 
% L8_20160604_pan.mat 
% L8_20160604_NDWI_SGL_mask.mat 

% Note that you already need a NDWI mask of lake locations (or indeed other
% method of getting lake areas if wished). Also needs a mask to show on-
% and off-ice sheet pixels (1 for on-ice, 0 for off-ice); reads in these at
% 10 m and 30 m resolution here. The 10 m resolution stuff can be commented
% out below if you don't want to resample the depth arrays to 10 m once
% lake depth calculations are complete. 

% Landsat-8 band files need to have dark offshore pixels visible in them
% for this to work, else you can set the R_inf values to 0 and it doesn't
% have much effect. 

% Andrew Williamson, University of Cambridge, August 2017 

%% Specify input file locations 

L8_main_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/'; 
L8_folder = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_LS8';
L8_bands_folder = [L8_folder '/bands'];
L8_NDWI_folder = [L8_folder '/NDWI_masks'];
L8_mask_folder = [L8_main_dir '/masks']; 
L8_out_depth_folder = [L8_folder '/depths']; 
L8_out_info_folder = [L8_folder '/lake_info']; 
L8_out_depth10_folder = [L8_folder '/depths_10m']; 
LS8_pix_siz = 0.03; LS8_pix_siz_m = 30; % specifying L8 resolution  
g_red = 0.7507; % value from Pope and others (2016) 
g_pan = 0.3817; % value from Pope and others (2016) 

%% Make directories  

% make new directories for the depth arrays (at original L8 30 m resolution
% and 10 m S2 resolution for analysis wrt S2) 
mkdir(L8_out_depth_folder); mkdir(L8_out_depth10_folder); 
mkdir(L8_out_info_folder); 

%% Load data into workspace 

cd(L8_bands_folder); 
L8_pan_list = dir('*_pan.mat'); 
L8_red_list = dir('*_red.mat'); 
cd(L8_NDWI_folder); 
L8_NDWI_list = dir('*_NDWI_SGL_mask.mat'); 
cd(L8_mask_folder); 
load('S2_mask_30m.mat'); 
GrIS_mask = S2_mask_30m; clear S2_mask_30m % variable here was saved as S2_mask_30m
GrIS_mask_10m = imread('S2_mask_10m.tif'); 

num_img = length(L8_pan_list); 
string = num2str(num_img); 
disp(['Calculating lake depths for ' string ' Landsat-8 images.'])

for i = 1:num_img 
    
    % get data and write image details to command line 
    dates = {L8_red_list(:).name}; 
    dates = dates'; % transpose matrix 
    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 = L8_red_list(i).name; 
    pan_fname = L8_pan_list(i).name; 
    NDWI_fname = L8_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(L8_bands_folder); 
    disp('Loading panchromatic band surface reflectance data.')
    load(pan_fname); 
    pan_refl = pan_band; 
    clear pan_band 
    disp('Loading red band surface reflectance data.')
    load(red_fname); 
    red_refl = red_band; 
    clear red_band 
    cd(L8_NDWI_folder); 
    disp('Loading NDWI-derived lake area mask.')
    load(NDWI_fname);
    lake_mask = NDWI_defined_SGLs; 
    %clear NDWI_defined_SGLs 
    
    % convert surface reflectance to double format for analysis and mask
    % no-data values (negatives in LS8) 
    pan_refl = double(pan_refl); 
    red_refl = double(red_refl); 
    pan_refl(pan_refl < 0) = NaN; 
    red_refl(red_refl < 0) = NaN; 
    % could potentially remove the above lines as they are now probably
    % done within the L8_preprocess script 
    
    % interpolate pan band to 30 m resolution (bilinearly following Pope
    % and others, 2016) 
    target_size = size(red_refl); 
    pan_refl = imresize(pan_refl,target_size,'bilinear'); 
    
    % get darkest pixel values for the images for calibrating the
    % depth-reflectance curve (see Williamson and others, 2017, for
    % details, based on Sneed and Hamilton, 2007)
    rinf_pan = min(min(pan_refl)); 
    rinf_red = min(min(red_refl)); 
    rinf_pan_str = num2str(rinf_pan,'%.3f'); 
    rinf_red_str = num2str(rinf_red,'%.3f'); 
    disp(['Darkest red band pixel value is ' rinf_red_str '.']); 
    disp(['Darkest panchromatic band value is ' rinf_pan_str '.']);
    
    if rinf_red > 0.1 
        disp('Warning: darkest red band value is high - manually check.'); 
    elseif rinf_pan > 0.1 
        disp('Warning: darkest panchromatic band value is high - manually check.'); 
    end
    
    pan_refl(GrIS_mask == 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; % remove NaNs by setting as 0s to 
    % allow the bwlabel to work 
    [lake_mask_labelled,lakenum] = bwlabel(lake_mask,8); % default 8-connectivity 
    
    % set up blank arrays for depth information 
    lake_depths = zeros(target_size); 
    lakeinfo = zeros(lakenum,9); 
    e = strel('disk',1); % to dilate lake mask by 1 pixel 
    
    % do the lake depth calculations for each lake in the array
    disp(['Calculating depths for all ' 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_pan = pan_refl; % tmpref is equal to values of pan band 
    tmpref_red = red_refl; % 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)) * LS8_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_pan = (pan_refl(ledge == 1)); % find the pan edge pixel refl
    edgepix_red = (red_refl(ledge == 1)); % find the red edge pixel refl 
    ad_pan = nanmean(edgepix_pan); % take the mean of these values (ignore NaNs) 
    ad_red = nanmean(edgepix_red); % take the mean of these values (ignore NaNs) 
    lakeinfo(n,5) = ad_red; % put this in lake info file 
    lakeinfo(n,9) = ad_pan; % put this in lake info file 
    
    % mask reflectance arrays 
    tmpref_pan(lmask == 0) = 0; 
    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_pan = (log(ad_pan - rinf_pan) - log(tmpref_pan - rinf_pan)) / (-g_pan);    
    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_pan; 
    tmpdepth = tmpdepth./2; % divide by two (i.e. take average)  
    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)) * (LS8_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  
    
    % make 10 m resolution depth array for Sentinel analysis too 
    S2_size = target_size * 3; 
    lake_depths_10m = imresize(lake_depths,S2_size,'nearest'); 
    
    % mask the lake depth arrays off the ice-sheet edge with NaNs 
    lake_depths_10m(GrIS_mask_10m == 0) = NaN; 
    lake_depths(GrIS_mask == 0) = NaN; 
    
    % get new variable names from above 
    new_depth_var_name = ['L8_' short_dt_str '_depths'];
    new_depth_10m_var_name = ['L8_' short_dt_str '_depths_10m'];
    new_info_var_name = ['L8_' short_dt_str '_lake_info'];
    
    % save the output files 
    cd(L8_out_depth_folder); 
    save([new_depth_var_name '.mat'],'lake_depths')
    cd(L8_out_depth10_folder); 
    save([new_depth_10m_var_name '.mat'],'lake_depths_10m')
    cd(L8_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 
    disp(['Lake depth calculations now ' num2str(progress,'%.1f') '% complete.']) 
    
end 

cd(L8_main_dir) 
