
% Script to set up maximum extent mask to show superimposed lake areas from
% each day of year within which a tracking algorithm can then be run. 

% Requires previous outputs of NDWI masks for each day of the year, and can
% be run either for Landsat-8 or Sentinel-2 individually or for both
% together; if run together and there are overlapping days with both image
% types available, then Sentinel-2 images are preferentially chosen. The
% NDWI masks must be the 10 m resolution ones.

% Andrew Williamson, October 2017 

%% Hard code parameters and directories 
 
clear all; clc; 
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';
tracking_out_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_lake_tracking_2017';
max_extent_dir = [tracking_out_dir '/max_extent/']; 
masks_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/masks'; 
radius = 200; % radius for eroding GrIS mask by (200 = 2 km) 
analysis_type = 'Sentinel'; % either 'Sentinel', 'Landsat' or 'both'
% depending on which images you'd like to use to make the mask from... 
duplicates = [20160926;20160731;20160701]; % dates with both Landsat-8 and Sentinel-2 images 
duplicates = []; % if there are no duplicates, so specify an empty array with [] 
mkdir(tracking_out_dir); mkdir(max_extent_dir); 

%% Load GrIS masking image & erode to remove erroneous lakes at ice margin

cd(masks_dir); 
GrIS_mask = imread('S2_mask_10m.tif'); 
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
% comment out two lines above if no erosion of mask required 

%% Set up blank maximum extent 

max_lake_extent = zeros(size(GrIS_mask)); 
% array of entirely 0s with the same size as the GrIS mask (which is the
% same size as the NDWI masks) 

%% Get list of images to add to mask 

cd(L8_lk_masks_dir);
L8_NDWI_list = dir('*_NDWI_10m_SGL_mask.mat'); 
cd(S2_lk_masks_dir); 
S2_NDWI_list = dir('*_NDWI_SGL_mask.mat'); 

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

%% Tidy up image list to remove duplicate dates if exist 

% Currently set up to deal with only 3 duplicate image dates, so the logic
% below will need altering if this changes 

% Currently also set up to remove the Landsat 8 imagery, not the Sentinel-2
% imagery, so if this changes then alter the below where 'L8' is specified 

if ~isempty(duplicates) 
    
img_num = length(images_to_add); 

dup_dt1 = duplicates(1,1); 
dup_dt2 = duplicates(2,1); 
dup_dt3 = duplicates(3,1); 

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,4:7);
    month = current_date(1,8:9); 
    day = current_date(1,10:11);
    image_type = current_date(1,1:2); 
    curr_date_str_shrt = [year month day]; 
    curr_date_num = str2double(curr_date_str_shrt); 
    date_nums(i,2) = curr_date_num; 
    date_nums(i,1) = i; 
    
    test1 = curr_date_num - dup_dt1; 
    test2 = curr_date_num - dup_dt2; 
    test3 = curr_date_num - dup_dt3; 
    
    if (test1 == 0 || test2 == 0 || test3 == 0) && strcmp(image_type,'L8') == 1 
        date_nums(i,3) = 1; 
    else 
        date_nums(i,3) = 0; 
    end 
    
end 

problems = date_nums; 
problems(problems(:,3) == 0,:) = []; 
problems = flipud(problems); 

for i = 1:length(problems) 
    
    row_num = problems(i,1); 
    images_to_add(row_num) = []; 
    
end 

end 

%% Add lakes from individual images to maximum extent mask 

% This is using the tidy list of imagery with any duplicates already removed 

img_num = length(images_to_add);

for i = 1:img_num 
    
    name = {images_to_add(i).folder}; 
    directory = name{1}; 
    date = {images_to_add(i).name}; 
    date = date'; 
    date = date{1}; 
    year = date(1,4:7); 
    month = date(1,8:9); 
    day = date(1,10:11); 
    image_type = 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'; 
    elseif strcmp(image_type,'S2') == 1 
        satellite = 'Sentinel-2'; 
    end 
    
    disp(['Adding lakes from ' curr_date_str ' from ' satellite ' to maximum lake extent mask.']) 
    pct_comp = 100/img_num * i; 
    disp(['Processing now ' num2str(pct_comp,'%.1f') '% complete.']); 
    
    cd(directory); 
    fname = images_to_add(i).name; 
    load(fname); 
    if strcmp(image_type,'L8') == 1 
        lakes = L8_lmask_10m; % if variable names changed in earlier lake 
        % area scripts, need to change here too, i.e. L8_lmask_10m 
    elseif strcmp(image_type,'S2') == 1 
        lakes = NDWI_defined_SGLs; 
    end 
    
    % remove cloud shadows from problem images for Landsat 8... 
    if curr_date_num == 20160716 && strcmp(image_type,'L8') == 1 
        lakes(4100:6500,9400:10600) = 0; 
        lakes(1:800,5600:6400) = 0;    
    elseif curr_date_num == 20160802 && strcmp(image_type,'L8') == 1 
        lakes(7800:8200,4800:5600) = 0; 
    end 
    
    % mask lakes with GrIS mask 
    lakes(GrIS_mask == 0) = 0; 
    
    % add these locations to the maximum extent mask
    max_lake_extent(lakes == 1) = 1;  
    
end 

%% View maximum lake extent figure 

figure; 
imagesc(max_lake_extent); 
title(['Maximum lake extent from ' analysis_type]); 

%% Save maximum extent mask for further use later 

cd(max_extent_dir); 
disp('Saving maximum lake extent.') 

if strcmp(analysis_type,'Sentinel') == 1 
    save('S2_max_lake_extent.mat','max_lake_extent'); 
elseif strcmp(analysis_type,'Landsat') == 1 
    save('L8_max_lake_extent.mat','max_lake_extent'); 
elseif strcmp(analysis_type,'both') == 1 
    save('S2_L8_max_lake_extent.mat','max_lake_extent');
end

cd(root_dir); 
