
% Script to pre-process Sentinel-2 data. 

% Sentinel-2 (10-m resolution) band specifications 
% Band 2 = blue
% Band 3 = green 
% Band 4 = red 

% S-2 image bands were pre-downloaded from the Sentinel Amazon Web Service
% and saved into the appt. directory, with bands having the following
% naming convention (this must be adhered to strictly for this script to 
% work): 
% YR-MO-DY-BNO.jp2, 
% where YR = year (e.g. 16), MO = month (e.g. 08), DY = day (e.g. 31) and
% NO = band number (e.g. 02). Example file name: 16-05-15-B03.jp2 

% N.B. No cloud masking capability within this script - must be done
% separately here 

% Andrew Williamson, Unviersity of Cambridge, August 2017 

clc; close all; clear all; 

%% Hard-code parameters and directories 

S2_main_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/'; 
S2_in_folder = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/raw_S2_images_2017';
S2_out_folder = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_S2_2017';
S2_out_bands_folder = [S2_out_folder '/bands'];
S2_out_RGBs_folder = [S2_out_folder '/RGBs'];
S2_out_RGB_images_folder = [S2_out_folder '/RGB_images'];
S2_out_NDWI_folder = [S2_out_folder '/NDWI_masks'];
S2_out_no_data_masks_folder = [S2_out_folder '/no_data_masks'];
NDWI_threshold = 0.25; % set the NDWI threshold for classifying SGLs 
% Should be defined based on sensitivity testing with S2_NDWI_sensitivity.m
% script 
save_to_wkspace = 'N'; % save variables to workspace? 'Y' or 'N' if
% specifying 'Y', if lots of files may crash... 
tidy_wkspace = 'Y'; % do you want to tidy up workspace at the end of the 
% script to just leave the variables set outside of for-loop? recommend 'Y'
% unless conducting testing 

%% Load Sentinel-2 data 

mkdir(S2_out_folder); mkdir(S2_out_bands_folder); mkdir(S2_out_RGBs_folder);
mkdir(S2_out_NDWI_folder); mkdir(S2_out_RGB_images_folder); 
mkdir(S2_out_no_data_masks_folder); 
cd(S2_in_folder); 
S2_b2_list = dir('*B02.jp2'); % blue band list
S2_b3_list = dir('*B03.jp2'); % green band list
S2_b4_list = dir('*B04.jp2'); % red band list
num_S2_files_b2 = length(S2_b2_list); % number of band 2 files 
num_S2_files_b3 = length(S2_b3_list); % number of band 3 files 
num_S2_files_b4 = length(S2_b4_list); % number of band 4 files 
string = num2str(num_S2_files_b2); % convert to string for command line 

disp(['In total, there are ' string ' Sentinel-2 images for analysis.']) 

check1 = num_S2_files_b2 - num_S2_files_b3; % checking number of b3 files 
check2 = num_S2_files_b2 - num_S2_files_b4; % checking number of b4 files 

% display status of checks at command line 
if check1 ~= 0 
    disp('Error! The number of band 3 files does not match the number of band 2 files.'); 
elseif check2 ~= 0 
    disp('Error! The number of band 4 files does not match the number of band 2 files.'); 
else
    disp('Good: correct number of files found for all bands. Continuing...');
end 
clear check1 check2 % no longer need these variables 
    
for i = 1:num_S2_files_b2
     
    % first get the date information from the input file
    % only need to do it for the band 2 since it should match for all
    % anyway, assuming they are in the file lists in the same order 
    dates = {S2_b2_list(:).name}; 
    dates = dates'; % transpose matrix 
    current_date = dates{i,1}; 
    year = current_date(1,1:2); %year_num = str2double(year); 
    month = current_date(1,4:5); %month_num = str2double(month); 
    day = current_date(1,7:8); %day_num = str2double(day); 
    curr_date_str = [day '/' month '/' year]; % format date string nicely 
    curr_imnum_str = num2str(i); 
    tot_im_str = num2str(num_S2_files_b2); 
    
    disp(['Now analysing image from ' curr_date_str ', image number '...
        curr_imnum_str ' of ' tot_im_str '.']); 
    
    % new variable for short name st
    short_dt_str = ['20' year month day]; 
    
    % find out which variable to read in
    blue_band_fname = S2_b2_list(i).name; 
    green_band_fname = S2_b3_list(i).name; 
    red_band_fname = S2_b4_list(i).name; 
    
    % read in the raw band files and assign variables 
    cd(S2_in_folder); 
    disp('Loading blue band data.') 
    blue_band = imread(blue_band_fname); 
    blue_band_forRGB = blue_band; 
    blue_band = double(blue_band); % convert to double for processing
    disp('Loading green band data.') 
    green_band = imread(green_band_fname); 
    green_band_forRGB = green_band; 
    green_band = double(green_band); % convert to double for processing
    disp('Loading red band data.') 
    red_band = imread(red_band_fname); 
    red_band_forRGB = red_band; 
    red_band = double(red_band); % convert to double for processing
        
    % make a mask to show the locations that have no data for them 
    no_data_mask = ones(size(red_band)); 
    no_data_mask(red_band == 0) = 0; 
    
    % show RGB image and stretch to the limits of the data 
    RGB_out = cat(3,red_band_forRGB,green_band_forRGB,blue_band_forRGB); 
    low_high = stretchlim(RGB_out); 
    new_RGB_image = imadjust(RGB_out,low_high,[]); 
    figure(1); 
    imshow(new_RGB_image); 
    cd(S2_out_RGB_images_folder);
    disp('Saving RGB image.') 
    out_im_name = ['S2_' short_dt_str '_RGB.jp2'];
    imwrite(new_RGB_image,out_im_name,'jp2','CompressionRatio',80,'Mode','lossy')
    
    % classify SGLs for this image with the NDWI technique 
    NDWI_str = num2str(NDWI_threshold); 
    disp(['Classifying SGLs based on the NDWI technique (threshold = ' NDWI_str ').'])  
    NDWI_defined_SGLs = S2_NDWI(blue_band,red_band,NDWI_threshold); % NDWI function 
    
    % filter small lakes (< 5 pixels) and linear lakes using function 
    NDWI_defined_SGLs = SGL_mask_filter(NDWI_defined_SGLs,5,2);
    
    % mask out no-data values from NDWI mask
    NDWI_defined_SGLs(no_data_mask == 0) = NaN; 
    
    % set new variable names with the short date string 
    new_blue_var_name = ['S2_' short_dt_str '_blue']; 
    new_green_var_name = ['S2_' short_dt_str '_green']; 
    new_red_var_name = ['S2_' short_dt_str '_red']; 
    new_RGB_var_name = ['S2_' short_dt_str '_RGB'];
    new_NDWI_var_name = ['S2_' short_dt_str '_NDWI_SGL_mask'];
    new_no_data_var_name = ['S2_' short_dt_str '_no_data_mask'];
    
    % save the bands to the output folders and remove from workspace 
    cd(S2_out_bands_folder);
    disp('Saving blue band.') 
    save([new_blue_var_name '.mat'],'blue_band')
    disp('Saving green band.') 
    save([new_green_var_name '.mat'],'green_band')
    disp('Saving red band.') 
    save([new_red_var_name '.mat'],'red_band')
    cd(S2_out_NDWI_folder); 
    disp('Saving NDWI mask.') 
    save([new_NDWI_var_name '.mat'],'NDWI_defined_SGLs')
    cd(S2_out_RGBs_folder); 
    disp('Saving RGB array.') 
    save([new_RGB_var_name '.mat'],'new_RGB_image')
    cd(S2_out_no_data_masks_folder);
    disp('Saving no-data mask.')
    save([new_no_data_var_name '.mat'],'no_data_mask'); 
    
    % rename variables appropriately with the date and the band number 
    % only if specified that you want to do this above 
    if strcmp(save_to_wkspace,'Y') == 1 
        
        assignin('base',new_blue_var_name,blue_band); 
        assignin('base',new_green_var_name,green_band);
        assignin('base',new_red_var_name,red_band);
        assignin('base',new_RGB_var_name,new_RGB_image);
        assignin('base',new_NDWI_var_name,NDWI_defined_SGLs);
        assignin('base',new_no_data_var_name,no_data_mask);
        
    end
    
    % get percentage progress update for command line and post to it 
    pct_complete = (100/num_S2_files_b2) * i; 
    pct_str = num2str(pct_complete); 
    disp(['Image processing is now ' pct_str '% complete.']); 
    disp(' '); 
    
    if strcmp(tidy_wkspace,'Y') == 1 
        
       clear blue_band red_band green_band new_blue_var_name ...
        new_red_var_name new_green_var_name red_band_fname ...
        green_band_fname blue_band_fname new_RGB_image RGB_out low_high ...
        curr_date_str current_date day day_num i month month_num ...
        new_RGB_var_name short_dt_str string year year_num ...
        NDWI_defined_SGLs new_NDWI_var_name blue_band_forRGB ...
        red_band_forRGB green_band_forRGB pct_str pct_complete ...
        no_data_mask out_im_name cmask_name curr_imnum_str ...
        new_no_data_var_name 
    % keeping workspace tidy (and small enough to allow processing!)
    
    end
    
end 

%% Do cloud masking for Sentinel-2 images 

% do the cloud masking for the images that have clouds present on them
% (which are specified as an input within this other script 
% full details of what this script does are provided within it
disp('Beginning to mask clouds from images with SWIR band.') 
run('S2_cloud_filter.m') 

%% Finish up 

cd(S2_main_dir); 
disp('All image processing complete.') 
