 
% Script to create cloud masks from Sentinel-2 imagery and then apply these
% to the NDWI masks to ensure that they are correctly filtered and no-data
% values are marked (instead of as 0s, which would cause misclassification
% of lake drainage further down the line). 

% Currently, the automatic cloud masking capabilities for Sentinel-2 are
% poor, meaning that some other method must be used. This is therefore a
% semi-automatic method. 

% Andrew Williamson, August 2017 

%% Hard code directories and parameters 

year_ref = '2017'; % change based on whether 2016 or 2017 analysis

S2_main_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/'; 
S2_out_folder = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_S2_2017';
S2_out_NDWI_folder = [S2_out_folder '/NDWI_masks'];
%S2_out_cmasked_NDWIs_folder = [S2_out_folder '/NDWI_masks_noclouds']; % folder for testing only 
cmasks_folder = [S2_out_folder '/cloud_masks']; 
S2_SWIR_bands = [S2_main_dir 'raw_S2_SWIR_2017']; 

if strcmp(year_ref,'2016') == 1 
    
img_checks = ['0502';'0601';'0604';'0608';'0611';'0617';'0618';...
    '0621';'0624';'0701';'0718';'0810';'0906']; 
% input list of the S-2 image dates that require checking for clouds (just
% use month and day here 

% these are included here because they are used to make arrays entirely of
% 0s as they are required to use the cloud mask arrays within the lake
% tracking approach, so i'm making a list of all of the images that are
% definitely cloud-free (100%) which are used below 
non_cloudy_imgs = ['0505';'0515';'0518';'0522';'0529';'0628';'0707';...
    '0708';'0711';'0714';'0721';'0724';'0727';'0728';'0731';'0803';...
    '0806';'0807';'0813';'0817';'0830';'0902';'0909';'0912';'0915';'0926'];

elseif strcmp(year_ref,'2017') == 1 
    
    img_checks = ['0501';'0503';'0504';'0510';'0513';'0520';'0521';...
        '0523';'0524';'0527';'0530';'0531';'0602';'0603';'0606';...
        '0610';'0612';'0613';'0616';'0619';'0626';'0629';'0701';'0703';...
        '0706';'0707';'0709';'0711';'0712';'0713';'0716';'0717';'0718';...
        '0721';'0722';'0723';'0724';'0727';'0728';'0729';'0730';'0731';...
        '0801';'0802';'0803';'0805';'0806';'0807';'0808';'0809';'0810';...
        '0811';'0812';'0813';'0814';'0815';'0820';'0821';'0822';'0824';...
        '0825';'0826';'0827';'0904';'0905';'0906';'0907';'0908';'0911';...
        '0912';'0918';'0920';'0921';'0922';'0923';'0924';'0925';'0926';...
        '0927';'0929';'0930']; 
    
    non_cloudy_imgs = [];   
end 

save_to_wkspace = 'N'; % save variables to workspace? 'Y' or 'N' if
% specifying 'Y', if lots of files may crash... 

%mkdir(S2_out_cmasked_NDWIs_folder); 
mkdir(cmasks_folder); 

%% Read in NDWI masks and SWIR bands for those images with clouds present 

% Not bothering with those images that are cloud free obviously... 

num_checks = length(img_checks); % find out how many images it's needed for
cd(S2_SWIR_bands); 

if strcmp(year_ref,'2016') == 1
    
    SWIR_list = dir('*-SWIR.jp2'); % get list of the SWIR bands (must have
    % this naming convention) 

elseif strcmp(year_ref,'2017') == 1 

    SWIR_list = dir('*-B11.jp2'); % get list of the SWIR bands (must have
    % this naming convention)  

end 

for i = 1:num_checks % for each image... 
    
    % go to the NDWI directory 
    cd(S2_out_NDWI_folder);  
    date = img_checks(i,:); % get the image date from the manual input
    mo = date(1:2); dy = date(3:4); num_str = num2str(num_checks); 
    curr = num2str(i); 
    disp(['Filtering clouds on ' dy '/' mo '/2016, image ' curr ...
        ' of ' num_str ' images requiring cloud masking.']); 
    
    if strcmp(year_ref,'2016') == 1 
        
    full_name = ['S2_2016' date '_NDWI_SGL_mask.mat']; 
    % specify name from this (based on the naming convention in the folder)
    
    elseif strcmp(year_ref,'2017') == 1 
        
    full_name = ['S2_2017' date '_NDWI_SGL_mask.mat']; 
    % specify name from this (based on the naming convention in the folder)
    
    end 
    var_name = ['lmask_' date]; % make new variable with this name 
    load(full_name); 
    
    if strcmp(save_to_wkspace,'Y') == 1 
        assignin('base',var_name,NDWI_defined_SGLs); % new variable assignment
    end
    
    cd(S2_SWIR_bands); % change to SWIR bands folder 
    SWIR_name = SWIR_list(i).name; % get the right name from the list above
    band11 = imread(SWIR_name); % read in image 
    band11 = double(band11); % convert to double for image processing 
    
    mo = SWIR_name(4:5); 
    dy = SWIR_name(7:8); 
    SWIR_var_name = ['SWIR_' mo dy]; % set up new variable name 
    
    tgt_sz = size(band11); 
    tgt_sz = tgt_sz .* 2; % SWIR bands are at 20 m resolution 
    
    band11 = imresize(band11,tgt_sz,'nearest'); % get to 10 m resolution 
    
    if strcmp(save_to_wkspace,'Y') == 1 
        assignin('base',SWIR_var_name,band11); % new variable 
    end
    % set up new array for cloud_mask with same size as band11
    cloud_mask = zeros(size(band11)); 
    cloud_mask(band11 > 1400) = 1; % mask out values where digital numbers
    % > 1400, the threshold based on manual comparison
    
    se = strel('disk',20);
    cloud_mask = imdilate(cloud_mask,se); 
    
    NDWI_defined_SGLs(cloud_mask == 1) = NaN; 
    
    cloudfree_lmask = NDWI_defined_SGLs; 
    
    cfree_var_name = ['lmask_cfree_' date];
    
    if strcmp(save_to_wkspace,'Y') == 1 
        assignin('base',cfree_var_name,cloudfree_lmask);
    end
    
    cd(S2_out_NDWI_folder) 
    save(full_name,'NDWI_defined_SGLs'); 
    
    cd(cmasks_folder);
    cmask_name = ['S2_2016' date '_cloud_mask.mat']; 
    save(cmask_name,'cloud_mask'); 
    
    pct_complete = (100/num_checks) * i; 
    pct_str = num2str(pct_complete); 
    disp(['Cloud masking is now ' pct_str '% complete.']); 
    disp(' '); 
    
end 

%% Save arrays of 0s for the images that are definitely cloud-free 

% This is to allow them to be used as part of the lake tracking process as
% it needs them to be able to run automatically and mask out clouds with
% these arrays from lake depth arrays etc. 

disp('Saving non-cloudy blank arrays for lake-tracking methods.') 

if ~isempty(non_cloudy_imgs) 
    
num_nc_img = length(non_cloudy_imgs); 

array_size = size(cloud_mask); 
cloud_mask = zeros(array_size); % nah it is used below (under save function)

for i = 1:num_nc_img 
    
    nc_date = non_cloudy_imgs(i,:); 
    
    cd(cmasks_folder);
    cmask_name = ['S2_2016' nc_date '_cloud_mask.mat']; 
    save(cmask_name,'cloud_mask'); 
    
end 

end 
%% Tidying up... 

clear NDWI_defined_SGLs full_name var_name date i num_checks ...
    band11 SWIR_var_name SWIR_name SWIR_list mo dy i tgt_sz cloud_mask se ...
    cloudfree_lmask cfree_var_name cmask dy mo pct_complete pct_str curr ...
    num_str 

cd(S2_main_dir); 
