
% Script to track lakes over series of Sentinel-2 MSI and Landsat 8 imagery
% using methods derived from Williamson et al.'s (2017) FAST algorithm, and
% quite similar to those used in Miles et al. (2017) with Sentinel-1 SAR and
% Landsat-8 OLI imagery. 

% Set up to allow running individually for Landsat-8, Sentinel-2 or both
% simultaneously, to see how outputs vary depending on this. 

% Andrew Williamson, University of Cambridge, 2017 

%% Hard-code parameters and directories 

root_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/'; 
L8_lk_masks_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_LS8/NDWI_10m_masks';
S2_lk_masks_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_S2_2017/NDWI_masks';
L8_lk_depth_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_LS8/depths_10m';
S2_lk_depth_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_S2_2017/depths'; 
L8_cld_mask_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_LS8/cloud_masks'; 
S2_cld_mask_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_S2_2017/cloud_masks'; 
L8_nan_mask_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_LS8/no_data_masks'; 
S2_nan_mask_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_S2_2017/no_data_masks'; 
tracking_out_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_lake_tracking_2017'; 
masks_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/masks'; 
max_extent_dir = [tracking_out_dir '/max_extent/']; 
imagery = 'Sentinel'; % can be set up for either Landsat-8, Sentinel-2 or
% both together: choose options 'Landsat', 'Sentinel' or 'both' 
analysis_yr = '2017'; % which year's data am I analysing? '2016' or '2017'
final_out_dir = [tracking_out_dir '/' imagery '/']; 
% make new directory for each type of imagery
radius = 200; % radius for eroding GrIS mask by (200 = 2 km) 
tracking_threshold = 495; % number of pixels used for lake tracking
pix_size = 0.01; % pixel size in km 
pix_size_m = 10; % pixel size in metres  
duplicates = [20160926;20160731;20160701]; 
% dates with contemporaneous L8 and S2 imagery to be removed from tracking,
% enter as column vector of the form defined above... 
max_poss_days = 153; % the maximum number of days in a season for which an
% image could ever be available (153 is from 1 May to 30 September) 
tidy_wkspace = 'Y'; % gets rid of unnecessary variables from workspace for
% tidiness on completion, enter 'Y' if you want to or anything else if not 
identify_drainage = 'Y'; % do you want to identify lake drainage events too
% (enter 'Y' if yes or anything else if not)  

%% Initialise 

disp('Lake tracking program beginning...') 
mkdir(tracking_out_dir); mkdir(final_out_dir); 
cd(masks_dir); 
disp('Loading ice mask and eroding edge by 2 km.') 
GrIS_mask = imread('S2_mask_10m.tif'); 
GrIS_mask_original = GrIS_mask; % making a copy of this for using as the
% uneroded version 
e = strel('disk',radius); % to get rid of the issues at the ice edge 
GrIS_mask = imerode(GrIS_mask,e); % erode mask by radius value
analysis_type = imagery; % setting to same value as above defined 'imagery'

%% Load file lists for analysis 

if strcmp(analysis_type,'Sentinel') == 1 
    
    disp('Conducting tracking for Sentinel-2 imagery alone.')
    
    % lake masks
    cd(S2_lk_masks_dir); 
    S2_NDWI_list = dir('*_NDWI_SGL_mask.mat');
    lake_masks = S2_NDWI_list; 
    
    % lake depths 
    cd(S2_lk_depth_dir); 
    S2_depth_list = dir('*_depths.mat'); 
    lake_depths_lst = S2_depth_list; 
    
    % maximum lake extent mask
    cd(max_extent_dir); 
    load('S2_max_lake_extent.mat'); 
    % the loaded MATLAB variable is called max_lake_extent 
    
    % cloud masks 
    cd(S2_cld_mask_dir); 
    S2_cloud_list = dir('*_cloud_mask.mat'); 
    cloud_masks = S2_cloud_list; 
    
    % no-data masks 
    % note Sentinel-2 data masks have 1s for data and 0s for no-data! 
    % Landsat 8 has the opposite (correction applied for consistency in
    % main for-loop below) 
    cd(S2_nan_mask_dir); 
    S2_no_data_list = dir('*_no_data_mask.mat'); 
    no_data_masks = S2_no_data_list; 
    
elseif strcmp(analysis_type,'Landsat') == 1 
    
    disp('Conducting tracking for Landsat 8 imagery alone.')
    
    % lake masks 
    cd(L8_lk_masks_dir);
    L8_NDWI_list = dir('*_NDWI_10m_SGL_mask.mat'); 
    lake_masks = L8_NDWI_list; 
    
    % lake depths 
    cd(L8_lk_depth_dir); 
    L8_depth_list = dir('*_depths_10m.mat');
    lake_depths_lst = L8_depth_list; 
    
    % maximum lake extent mask
    cd(max_extent_dir); 
    load('L8_max_lake_extent.mat'); 
    
    % cloud masks 
    cd(L8_cld_mask_dir); 
    L8_cloud_list = dir('*_cloud_mask.mat'); 
    cloud_masks = L8_cloud_list; 
    
    % no-data masks 
    % note Sentinel-2 data masks have 1s for data and 0s for no-data! 
    % Landsat 8 has the opposite (correction applied for consistency in
    % main for-loop below) 
    cd(L8_nan_mask_dir); 
    L8_no_data_list = dir('*_no_data_mask.mat'); 
    no_data_masks = L8_no_data_list; 
    
elseif strcmp(analysis_type,'both') == 1
    
    disp('Conducting tracking for both Sentinel-2 and Landsat 8 imagery.')
    
    % lake masks 
    cd(S2_lk_masks_dir); 
    S2_NDWI_list = dir('*_NDWI_SGL_mask.mat');
    cd(L8_lk_masks_dir);
    L8_NDWI_list = dir('*_NDWI_10m_SGL_mask.mat'); 
    lake_masks = [S2_NDWI_list;L8_NDWI_list];
    
    % lake depths 
    cd(L8_lk_depth_dir); 
    L8_depth_list = dir('*_depths_10m.mat');
    cd(S2_lk_depth_dir); 
    S2_depth_list = dir('*_depths.mat'); 
    lake_depths_lst = [S2_depth_list;L8_depth_list]; 
    
    % maximum lake extent mask 
    cd(max_extent_dir); 
    load('S2_L8_max_lake_extent.mat'); 
    
    % cloud masks 
    cd(S2_cld_mask_dir); 
    S2_cloud_list = dir('*_cloud_mask.mat'); 
    cd(L8_cld_mask_dir); 
    L8_cloud_list = dir('*_cloud_mask.mat'); 
    cloud_masks = [S2_cloud_list;L8_cloud_list];
    
    % no-data masks 
    % note Sentinel-2 data masks have 1s for data and 0s for no-data! 
    % Landsat 8 has the opposite (correction applied for consistency in
    % main for-loop below) 
    cd(L8_nan_mask_dir); 
    L8_no_data_list = dir('*_no_data_mask.mat'); 
    cd(S2_nan_mask_dir); 
    S2_no_data_list = dir('*_no_data_mask.mat'); 
    no_data_masks = [S2_no_data_list;L8_no_data_list];
    
end 

%% Tidy up file lists to remove duplicate dates for Sentinel-2 & Landsat 8

% get the number of images from the length of the lake_masks list 
img_num = length(lake_masks); 

% extract the duplicate dates from the previously provided list (note
% column vector) 
dup_dt1 = duplicates(1,1); % first date 
dup_dt2 = duplicates(2,1); % second date 
dup_dt3 = duplicates(3,1); % third date 

% set up a blank array with same number of rows as img_num and 3 columns 
date_nums = zeros(img_num,3); 
dates = {lake_masks(:).name}; % get the names of the images from lake_masks
dates = dates'; % make vertical 

% for each of the dates in list, do this... 
for i = 1:img_num 
    
    current_date = dates{i,1}; 
    year = current_date(1,4:7);
    month = current_date(1,8:9); 
    day = current_date(1,10:11);
    image_type = current_date(1,1:2); % find out if "L8" or "S2" 
    curr_date_str_shrt = [year month day]; 
    curr_date_num = str2double(curr_date_str_shrt); % convert to double 
    date_nums(i,2) = curr_date_num; % put the date into 2nd column 
    date_nums(i,1) = i; % put the number into 1st column
    
    % do tests to check whether the current date is the same as the
    % duplicate dates... (3 times for 3 dates that overlap in 2016) 
    test1 = curr_date_num - dup_dt1; 
    test2 = curr_date_num - dup_dt2; 
    test3 = curr_date_num - dup_dt3; 
    
    % if the tests = 0 (i.e. image dates are same) and the imagery type is
    % Landsat-8 (the type we want to remove), then put a 1 in the 3rd
    % column, else put a 0 in it (i.e. don't remove it) 
    if (test1 == 0 || test2 == 0 || test3 == 0) && ...
            strcmp(image_type,'L8') == 1 
        date_nums(i,3) = 1; 
        
        % this is to override this test when doing only one type of imagery
        % to avoid screwing up the missing info required in the date_nums
        % array (needed later for a multitude of reasons) 
        if strcmp(imagery,'Landsat') == 1 || strcmp(imagery,'Sentinel') == 1
            date_nums(i,3) = 0; 
        end 
        
    else 
        date_nums(i,3) = 0; 
    end 
    
end 

% find out the problematic (i.e. duplicate) days 
problems = date_nums; % make a copy of this array 
problems(problems(:,3) == 0,:) = []; % get rid of those with 0s here as 
% we don't want to remove them 
problems = flipud(problems); % so that it removes the latest one first or
% else row references are incorrect 

% for each of the duplicated dates... 
for i = 1:length(problems) 
    
    % extract the row number of the duplicate  
    row_num = problems(i,1); 
    lake_masks(row_num) = []; % remove the row from the list 
    no_data_masks(row_num) = []; % same for the no-data masks list 
    lake_depths_lst(row_num) = []; % and lake depths list 
    date_nums(row_num,:) = []; % and date nums list 
    cloud_masks(row_num) = []; % and cloud masks list 
    
end 

img_num = length(lake_masks); 
temp_res = max_poss_days / img_num; 
disp(['In total, there are ' num2str(img_num) ' images for lake tracking.'])
disp(['This represents an average temporal image resolution of '...
    num2str(temp_res,'%1.1f') ' days over the entire sensing period.']) 

%% Label maximum extent 'regions' defined from all images 

tracking_mask = max_lake_extent; 
% set tracking mask same as max lake extent

labelled_track_mask = bwlabel(tracking_mask); 
% create new array and label it  
filledareas = regionprops(labelled_track_mask,'FilledArea'); 
% use regionprops functioned to find filledareas within this array 
filledareas2 = cat(1,filledareas.FilledArea); % concatenate 
lakes_to_track = find(filledareas2 >= tracking_threshold); 
% defines the lakes to track based on threshold set earlier in script. a
% connected region must be at least equal to or greater than this value to
% be examined over multiple images 

% creates filled lake mask and removes those elements that are not being
% tracked 
binary_image = ismember(labelled_track_mask,lakes_to_track); 
lake_mask_filled = labelled_track_mask; 
lake_mask_filled(~binary_image) = 0; % sets elements not part of image to 0 

% create labelled lake mask based upon the filled lake mask 
% using 8 connectivity
% tottracknum is also created, which is the total number of lakes to
% examine from maximum extent 
[lakemask_lab,tracking_num] = bwlabel(lake_mask_filled,8); 

disp(['Identified ' num2str(tracking_num)...
    ' lakes to track from maximum extent of all dates.']); 

% make figure of this, which has lakes colour-coded by their number 
figure; clf 
imagesc(lakemask_lab); 
title('Colour-coded maximum extent of lakes in 2016') 
hold on
contour(GrIS_mask_original,10,'k'); 

% save the labelled mask in case needed later 
cd(final_out_dir); 
save('lakes_tracked_labelled.mat','lakemask_lab'); 
save('lakes_inc_below_tracking_threshold.mat','labelled_track_mask'); 

%% Get boundaries of maximum extent 'regions' and make figure of it 

% don't need this section for analysis 

[Bo,Lo] = bwboundaries(lake_mask_filled,8,'holes'); 
% find boundaries based on holes within the image (again, 8-connectivity) 
j = 1;
lakes = 1;
boundaries = struct('lakes',lakes); % creates new structure 

% repeats this for-loop for each area within the image 
for j = 1:length(Bo)
    boundary = Bo{j};
    boundaries(j).lakes = boundary;
end

figure; clf;  
imagesc(lake_mask_filled); 
hold on 
contour(GrIS_mask_original,10,'k'); 

for n = 2:length(boundaries) % repeat for all boundaries 
    
    plot(boundaries(n).lakes(:,2),boundaries(n).lakes(:,1),'r',...
        'LineWidth', 1)
    % plot boundary for maximum extent as red lines with width 1 
    
    hold on % specify here to ensure all boundaries added sequentially 
    
end 

hold on

%% Find out maximum extent indices of lakes 

max_extent_indices = zeros(tracking_num,3); 

for n = 1:tracking_num % for the entire number of lakes 
           
    [col,row] = find(lakemask_lab == n); % x,y co-ordinates of lake
           
    % lines immediately below are for figure plotting purposes 
    textbox_x = max(row); % the x-coord to plot the left of the box
    textbox_y = mean(col); % same for y-coord
    text(textbox_x,textbox_y,sprintf(' %d',n),'FontSize',9,'Color','w');
    
    col_mean = mean(col); 
    row_mean = mean(row); 
    
    max_extent_indices(n,1) = n; 
    max_extent_indices(n,2) = row_mean; 
    max_extent_indices(n,3) = col_mean; 
    % Think these are the right way around in terms of row, col format! 
    
end

title('Labelled maximum extent of lakes in 2016 with boundaries') 

%% Calculate area and volume within each maximum extent region for each image

% ultimately this section results in an output of an x by y array, 
% where x is number of days of images and y is number of lakes examined 

area_info_file = zeros(size(tracking_num,img_num)); 
% create an array of zeros with the number of rows equal to the total
% number of lakes being tracked and the number of columns equal to the
% number of images for the year 

pixel_info_file = zeros(size(tracking_num,img_num)); 
% do the same, but this time for the number of pixels insead of area
% calculations 

volume_info_file = zeros(size(tracking_num,img_num)); 
% do the same, but this time for the volumes 

cloud_info_file = zeros(img_num,2); % and clouds
data_cover_info_file = zeros(img_num,2); % and data cover 

% for each image from the first to the total number part of the workspace,
% do the following: 
for nla = 1:img_num % img_num is total number of images 
    
    curr_area_array_name = lake_masks(nla).name;
    curr_area_array_dir = lake_masks(nla).folder; 
    curr_no_data_array_name = no_data_masks(nla).name; 
    curr_no_data_array_dir = no_data_masks(nla).folder; 
    curr_depth_array_name = lake_depths_lst(nla).name; 
    curr_depth_array_dir = lake_depths_lst(nla).folder; 
    curr_cloud_array_dir = cloud_masks(nla).folder; 
    curr_cloud_array_name = cloud_masks(nla).name; 
    image_type = curr_no_data_array_name(1,1:2); % find out if L8 or S2 
    cd(curr_area_array_dir); 
    load(curr_area_array_name); 
    cd(curr_no_data_array_dir); 
    load(curr_no_data_array_name); 
    cd(curr_depth_array_dir); 
    load(curr_depth_array_name); 
    cd(curr_cloud_array_dir); 
    load(curr_cloud_array_name); % has variable name of cloud_mask 
    
    if strcmp(image_type,'L8') == 1 
        curr_area_array = L8_lmask_10m; 
        curr_vol_array = lake_depths_10m; 
        no_data_mask = no_data_mask_10m; 
        no_data_mask(no_data_mask == 1) = 2; 
        no_data_mask(no_data_mask == 0) = 1;
        no_data_mask(no_data_mask == 2) = 0;
        tgt_siz = size(curr_vol_array); 
        cloud_mask = imresize(cloud_mask,tgt_siz,'nearest'); 
        % tgt_size based on correct resolution of curr_vol_array 
    elseif strcmp(image_type,'S2') == 1 
        curr_area_array = NDWI_defined_SGLs; 
        curr_vol_array = lake_depths; 
        % no_data_mask renaming not needed here as already called that 
        % same for cloud_mask 
    end  
    
    % this is to set the lake depth array to 0s wherever it currently has
    % NaN values (from the lake depth script), becuase otherwise this
    % messes with the calculations for lake volume tracking if there are
    % NaN values in it 
    logical_nan_vols = isnan(curr_vol_array); 
    curr_vol_array(logical_nan_vols == 1) = 0; 
    
    % no-data masks have 0s for no data, whereas cloud masks have 1s for no
    % data - super confusing... 
    curr_area_array(no_data_mask == 0) = NaN; 
    curr_vol_array(no_data_mask == 0) = NaN; 
    curr_area_array(cloud_mask == 1) = NaN; 
    curr_vol_array(cloud_mask == 1) = NaN; 
    
    % find out cloud cover and no data cover percentage relative to whole
    % region 
    dims = size(curr_area_array); 
    dims = dims(1,1) * dims(1,2); 
    dims = dims * pix_size * pix_size; % in km2 
    sum_clouds = sum(sum(cloud_mask)); 
    sum_no_data = sum(sum(no_data_mask)); 
    sum_clouds = sum_clouds * pix_size * pix_size; 
    sum_no_data = sum_no_data * pix_size * pix_size; 
    pct_cloud_cover = 100 / dims * sum_clouds; 
    pct_data_cover = 100 / dims * sum_no_data;
    cloud_info_file(nla,1) = pct_cloud_cover; 
    data_cover_info_file(nla,1) = pct_data_cover; 
    
    disp(['Calculating lake areas and volumes for image number ' ...
        num2str(nla) ' of ' num2str(img_num) ' total images.']) 
    
    % for each lake within this image from 1 until the total number
    % identified as part of the maximum extent array, do the following: 
    for nlg = 1:tracking_num 
        
        disp(['Examining lake number ' num2str(nlg) ' of ' ...
            num2str(tracking_num) ' lakes on the maximum extent.']); 
        % progress updates to check it's working as it takes a long time 
        
        curr_image_mask = zeros(size(curr_area_array)); 
        % make a new mask array of all zeros with same size as
        % curr_area_array 
        
        curr_image_mask(lakemask_lab == nlg) = 1; 
        % where the lakemask_lab array has a value equal to the nlg (i.e.
        % the current loop number), set the current image mask equal to 1; 
        % this means all areas outside the current lake are 0s and the
        % areas within the lake are 1s 
        
        curr_lake = curr_area_array(curr_image_mask == 1); 
        % this pulls out the current lake from the big mask by removing the
        % area from the current area array that is equal to 1 in the
        % current image mask 
        
        curr_lake_for_vol = curr_vol_array(curr_image_mask == 1); 
        % should technically be called "curr_depth_array" as these are
        % depths at this point, not volumes! can only calculate volumes
        % from the depths 
        
        curr_lake_pix = sum(sum(curr_lake)); 
        % find out the number of pixels in the lake area by summing the
        % number of elements that are exactly equal to 1 (i.e. == 1) in the
        % current lake 
        
        pixel_info_file(nlg,nla) = curr_lake_pix; 
        % write out to the pixel infor file for the correct day & correct
        % lake the number of pixels within that lake 
        
        curr_lake_area = curr_lake_pix * (pix_size ^ 2); 
        % work out area for the lake pixels by multiplying by the pixel
        % size squared 
        
        area_info_file(nlg,nla) = curr_lake_area; 
        % put this information into the lake area array
        
        curr_lake_vol = sum(sum(curr_lake_for_vol)); 
        % sum up all the depths to get 
        
        curr_lake_vol = curr_lake_vol * (pix_size_m ^ 2); 
        % getting lake volume from lake depth by multiplying by pixel size
        % squared 
        
        volume_info_file(nlg,nla) = curr_lake_vol; 
        
    end % end of nested for-loop  
    
    pct_comp = 100/img_num * nla; 
    disp(['Lake tracking process now ' num2str(pct_comp,'%1.1f') '% complete.']); 
    
end % end of main for-loop 

%% Remove infinite values from the arrays if present 

test_inf_vol = isinf(volume_info_file); 
volume_info_file(test_inf_vol == 1) = NaN; 
test_inf_area = isinf(area_info_file); 
area_info_file(test_inf_area == 1) = NaN; 
test_inf_pix = isinf(pixel_info_file); 
pixel_info_file(test_inf_pix == 1) = NaN; 

% because the lake depths are negative in 2017 (error with lake depth code,
% but it has no implications, simply a badly specified hard-coded
% variable... depths are all correct) 
if strcmp(analysis_yr,'2017') == 1 
    volume_info_file = -volume_info_file; 
end 

%% Derive array for total lake area and volume change over all imagery 

total_lake_areas = zeros(size(2,img_num)); 
% new array for total lake area for each image with 2 rows and the number
% of columns set according to the total number of images 

total_lake_volumes = zeros(size(2,img_num)); 

disp('Calculating total lake areas over season.') 

for i = 1:img_num % from 1 until the total number of images, do...
    
    curr_area_col = area_info_file(:,i); % extract the column from the area
    % info file as defined by the loop index and extract every row (:) 
    curr_vol_col = volume_info_file(:,i); 
    
    curr_area_date = date_nums(i,2); % extract date from the date_nums 
    % double file 
    curr_vol_date = date_nums(i,2); 
    
    total_lake_areas(1,i) = nansum(curr_area_col); 
    % put the sum of the lake areas into the first row of the
    % total_lake_areas array and the column to be defined the loop index 
    % nansum is used to take the sum ignoring any NaN values (in case this
    % is used earlier) 
    total_lake_volumes(1,i) = nansum(curr_vol_col); 
    
    total_lake_areas(2,i) = curr_area_date; % current JD number in 2nd row 
    total_lake_volumes(2,i) = curr_vol_date; 
end

total_lake_areas = total_lake_areas'; 
total_lake_volumes = total_lake_volumes'; 
total_lake_areas = sortrows(total_lake_areas,2); 
total_lake_volumes = sortrows(total_lake_volumes,2); 

total_cloud_cover = cloud_info_file; 
total_data_cover = data_cover_info_file; 
dates_only = date_nums(:,2); % get the date numbers only
dates_only = dates_only'; % make into a column 
total_cloud_cover(:,2) = dates_only; 
total_data_cover(:,2) = dates_only; 
total_data_cover = sortrows(total_data_cover,2); 
total_cloud_cover = sortrows(total_cloud_cover,2); 

cd(final_out_dir); % change to outputs directory 
save('total_lake_area.mat','total_lake_areas'); % save the lake area file 
save('total_lake_volume.mat','total_lake_volumes'); % volumes
save('total_cloud_cover.mat','total_cloud_cover'); % cloud cover
save('total_data_cover.mat','total_data_cover'); % data cover 

%% Add dates to pixel, area and volume info files for easier analysis 

dates_only = date_nums(:,2); % get the date numbers only
dates_only = dates_only'; % make into a row 

row_allocation = size(area_info_file,1)+1; % find out where to put the new info

% put the dates into the right place for the area array (then do the same
% for the volume and pixel arrays below too) 
area_info_file(row_allocation,:) = dates_only; 
area_info_file = area_info_file'; 
area_info_file = sortrows(area_info_file,row_allocation); 

volume_info_file(row_allocation,:) = dates_only; 
volume_info_file = volume_info_file'; 
volume_info_file = sortrows(volume_info_file,row_allocation); 

pixel_info_file(row_allocation,:) = dates_only; 
pixel_info_file = pixel_info_file'; 
pixel_info_file = sortrows(pixel_info_file,row_allocation); 

% swap rows to columns for all arrays 
area_info_file = area_info_file'; 
volume_info_file = volume_info_file'; 
pixel_info_file = pixel_info_file'; 

%% Finish 

cd(final_out_dir); % change to outputs directory 
save('lake_area_tracking.mat','area_info_file'); % save the lake area file
save('lake_pixel_tracking.mat','pixel_info_file'); % save the lake area file
save('lake_volume_tracking.mat','volume_info_file'); % save volume info file
cd(root_dir); % change back to main working directory 

disp('Lake tracking process done... Have fun!') 
cd(root_dir); 

if strcmp(identify_drainage,'Y') == 1 % if  wanting to identify drainage
    cd(root_dir); 
    run('S2_L8_lake_drainage.m'); % now run the lake tracking script
end 

% if wanting to tidy up workspace, then clear the following variables 
if strcmp(tidy_wkspace,'Y') == 1

    clear binary_image Bo boundaries boundary cloud_mask cloud_masks col ...
    curr_area_array curr_area_array_dir curr_area_array_name ...
    curr_cloud_array_dir curr_cloud_array_name curr_date_num ...
    curr_date_str_shrt curr_depth_array_dir curr_depth_array_name ...
    curr_image_mask curr_lake curr_lake_area curr_lake_for_vol ...
    curr_lake_pix curr_lake_vol curr_no_data_array_dir curr_no_data_array_name ...
    curr_vol_array current_date day dims dup_dt1 dup_dt2 dup_dt3 duplicates ...
    e filledareas filledareas2 i image_type j L8_cld_mask_dir L8_cloud_list ...
    L8_depth_list L8_lk_depth_dir L8_nan_mask_dir L8_NDWI_list L8_no_data_list ...
    lake_depths lake_depths_lst lake_masks lakes lakes_to_track Lo masks_dir...
    max_extent_dir max_lake_extent max_poss_days month n NDWI_defined_SGLs ...
    nla nlg no_data_mask no_data_masks pct_cloud_cover pct_data_cover problems ...
    radius row row_num S2_cld_mask_dir S2_cloud_list S2_depth_list S2_lk_depth_dir ...
    S2_nan_mask_dir S2_NDWI_list S2_no_data_list sum_clouds sum_no_data ...
    temp_res test1 test2 test3 textbox_x textbox_y tracking_mask tracking_num ...
    tracking_out_dir tracking_threshold year 
end 
