% (c) Rebecca Dell, University of Cambridge, June 2020.
% Adapted from the work of Williamson et al. (2018) (doi:10.5194/tc-2018-56). 
    
% Script to set up maximum extent mask to show superimposed water body areas from
% each day analysed within which a tracking algorithm can then be run. 

% Can run for either the 'combined' mask, the 'lake' masks or the 'stream'
% masks. 
% Can be run either for Landsat-8 or Sentinel-2 individually or for both
% together. The NDWI masks must be the 10 m resolution ones when running for both.
%% Hard code parameters and directories 
 
clear all; clc; 
root_dir = 'Dropbox (Cambridge University)'; 
L8_lk_masks_dir ='C:\Users\rld46\Dropbox (Cambridge University)\testing_code\CategorisedL8';
S2_lk_masks_dir = 'C:\Users\rld46\Dropbox (Cambridge University)\testing_code\CategorisedS2';
max_extent_dir = 'C:\Users\rld46\Dropbox (Cambridge University)\testing_code\tracking';
L8_mask_folder = 'C:\Users\rld46\Dropbox (Cambridge University)\testing_code\input_L8';
analysis_type = 'both'; % either 'Sentinel', 'Landsat' or 'both'

%% Get list of images to add to mask 

%NB. At the moment the script is designed to calculate the max extent mask
%for only one category at a time (e.g. combined, lakes, streams). This
%means variables have to be altered depending on what you want. 

cd(L8_lk_masks_dir);
L8_NDWI_combined_list = dir('*area_c.mat');
% L8_NDWI_lakes_list = dir('*_lakes_matrix.mat'); 
% L8_NDWI_streams_list = dir('*_streams_matrix.mat');
cd(S2_lk_masks_dir); 
S2_NDWI_combined_list = dir('*area_c.mat');
% S2_NDWI_lakes_list = dir('*_lakes_matrix.mat');
% S2_NDWI_streams_list = dir('*_streams_matrix.mat');

if strcmp(analysis_type,'Sentinel') == 1 
    images_to_add = S2_NDWI_list; 
elseif strcmp(analysis_type,'Landsat') == 1 
    images_to_add = [L8_NDWI_combined_list]
elseif strcmp(analysis_type,'both') == 1
    images_to_add = [L8_NDWI_combined_list;S2_NDWI_combined_list]; 
end 

%% Set up blank maximum extent 

cd(L8_mask_folder); 
load('ice_ocean_mask_S2.mat');
max_combined_extents = zeros(size(ice_ocean_mask_S2));

% Add lakes from individual images to maximum extent mask 

img_num = length(images_to_add); 
date_nums = zeros(img_num,3); 
dates = {images_to_add(:).name}; 
dates = dates'; 

for i = 1:img_num 

current_date = dates{i,1}; 
year = current_date(1,3:4);
month = current_date(1,5:6); 
day = current_date(1,7:8);
image_type = current_date(1,1:2);

    curr_date_str = [day '/' month '/' year]; % format date string nicely
    curr_date_str_shrt = [day month year]; 
    curr_date_num_pre = [year month day]; 
    curr_date_num = str2double(curr_date_num_pre); 
    if strcmp(image_type,'L8') == 1 
        satellite = 'Landsat 8'; 
          cd(L8_lk_masks_dir);
    elseif strcmp(image_type,'S2') == 1 
        satellite = 'Sentinel-2'; 
          cd(S2_lk_masks_dir);
    end 
    
    disp(['Adding water bodies from ' curr_date_str ' from ' satellite ' to maximum combined extents mask.']) 
    pct_comp = 100/img_num * i; 
    disp(['Processing now ' num2str(pct_comp,'%.1f') '% complete.']); 
    
    % Load the necessary mask
    fname = images_to_add(i).name; 
    load(fname); 
    
    % Add lakes and streams to the maximum extent mask
    %add these locations to the maximum extent mask
    if strcmp(image_type,'L8') == 1 
    max_combined_extents(curr_combined_area_array > 0)=1;
    elseif strcmp(image_type,'S2') == 1 
    max_combined_extents(curr_combined_area_array > 0)=1; 
    end
    
end
    
% View maximum lake extent figure 

figure; 
imagesc(max_combined_extents);colorbar; 
title(['Maximum combined extent from ' analysis_type]); 

%% Save maximum extent mask for further use later 

cd(max_extent_dir); 
disp('Saving maximum combined extents.') 

if strcmp(analysis_type,'Sentinel') == 1 
    save('S2_max_combined_extent.mat','max_combined_extents'); 
elseif strcmp(analysis_type,'Landsat') == 1 
    save('L8_max_combined_extent.mat','max_combined_extents'); 
elseif strcmp(analysis_type,'both') == 1 
    save('S2_L8_max_combined_extent.mat','max_combined_extents');
end







