
% Script to test Sentinel-2 sensitivity to NDWI thresholds by overlaying
% NDWI-delineated lake outlines onto background RGB images to allow the two
% to be compared and so to determine the best threshold. 
% For S-2, band 2 is blue and band 4 is red. 

% Andrew Williamson, University of Cambridge, July 2017 

%% Define input parameters for script 

S2_main_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/'; % root directory 
S2_in_folder = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/raw_images'; % raw S-2 bands
S2_out_folder = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs'; % folder for outputs 
S2_out_NDWI_sens_folder = [S2_out_folder '/NDWI_sensitivity']; % new directory to make in the script
S2_out_RGB_images_folder = [S2_out_folder '/RGB_images']; % previously created RGB images from S-2 bands
mask_folder = [S2_main_dir '/mask']; 
S2_mask_name = 'S2_mask_10m.tif'; 

NDWI_threshold = 0.36; % set the NDWI threshold for classifying SGLs 
img_dt = '2016-06-21'; % which image date are you doing this for? in YYYY-MO-DY format 
bands_loaded = 'N'; % cut out the loading if bands already in workspace; enter 'Y' or 'N' 

%% Start of main script 

if strcmp(bands_loaded,'N') == 1 % only do the below commands if bands not already loaded 
    
%% Load the ice-sheet mask, the blue and red bands and the RGB image 
% make directory if it doesn't already exist 
mkdir(S2_out_NDWI_sens_folder); 

% get the characters for reading file names 
yr_short = img_dt(3:4); yr = img_dt(1:4); 
mo = img_dt(6:7); dy = img_dt(9:10);

% load the ice sheet mask 
cd(mask_folder); 
mask = imread(S2_mask_name); 
disp('Loaded mask.') 

cd(S2_in_folder); % change to the input file location 
blue_str = [yr_short '-' mo '-' dy '-B02.jp2']; % make correct filename format for blue band 
blue_band = imread(blue_str); % read blue band
blue_band = double(blue_band); 
disp('Loaded blue band.') 
red_str = [yr_short '-' mo '-' dy '-B04.jp2']; % make correct filename format for red band 
red_band = imread(red_str); % read red band 
red_band = double(red_band); 
disp('Loaded red band.') 

cd(S2_out_RGB_images_folder); 
RGB_str = ['S2_' yr mo dy '_RGB.jp2']; % make correct filename format for RGB image 
RGB_im = imread(RGB_str); % read RGB image 
disp('Loaded RGB image.') 

%% Mask ice-marginal areas 

blue_band(mask == 0) = NaN; % set the blue band to NaN where the mask == 0
red_band(mask == 0) = NaN; 

end 

%% Do the NDWI testing 

NDWI_defined_SGLs = S2_NDWI(blue_band,red_band,NDWI_threshold); 
% a pre-existing function to keep the NDWI calculations separate 
% output from the above is the NDWI_defined_SGLs array, which will now be
% part of workspace for further analysis; this is a binary mask of lake
% locations (1s = lakes, 0s = not lakes) 
NDWI_thresh_str = num2str(NDWI_threshold); 

disp(['Calculated lake areas using the NDWI with threshold of ' NDWI_thresh_str]) 

%% Get lake boundaries from the NDWI 

lake_bounds = bwboundaries(NDWI_defined_SGLs,8,'holes'); % bwboundaries in-built function 
k = 1; lakes = 1; 
boundaries = struct('lakes',lakes); % make a structure containing each lake 

% get the boundary information 
for k = 1:length(lake_bounds) 
    boundary = lake_bounds{k}; 
    boundaries(k).lakes = boundary;
end

% get the lake mask and number of lakes in case need later 
[lakemask,lakenum] = bwlabel(NDWI_defined_SGLs,8); 

%% Calculate total lake area 

sum_pix = sum(sum(NDWI_defined_SGLs)); 
tot_lake_area = sum_pix * 0.01 * 0.01; 
tot_lk_area_str = num2str(tot_lake_area); 
disp(['Total lake area is ' tot_lk_area_str ' km^2']) 

%% Fig 1. The NDWI mask with boundaries overlaid to check it works 

figure(1) 
imagesc(NDWI_defined_SGLs); 
hold on 

% plot boundaries around the binary polygons of lakes 
for n = 1:length(boundaries) 
    plot(boundaries(n).lakes(:,2),boundaries(n).lakes(:,1),'r','LineWidth',1,'LineStyle','-')
    hold on 
end 

%% Fig 2. The RGB image with boundaries overlaid for manual comparison

figH = figure(2); 
set(figH,'name',['NDWI threshold of ' NDWI_thresh_str]); 
set(figH,'PaperType','a4','PaperOrientation','portrait','PaperUnits','normalized'); 
imshow(RGB_im) 
hold on 

% plot boundaries onto the RGB image 
for n = 1:length(boundaries) 
    plot(boundaries(n).lakes(:,2),boundaries(n).lakes(:,1),'r','LineWidth',1,'LineStyle','-')
    hold on 
end 
 
image_sens_folder = [S2_out_NDWI_sens_folder '/' img_dt];
mkdir(image_sens_folder); 
cd(image_sens_folder) 
figname = ['lakes_NDWI_' NDWI_thresh_str '.fig']; 
PDFname = ['lakes_NDWI_' NDWI_thresh_str '.pdf']; 
%saveas(gcf,PDFname,'pdf') 
%savefig(figname); 
disp(['Saved MATLAB figure window to ' image_sens_folder]); 

%% Fig 3. The RGB image with translucent array of NDWI lakes 

% This doesn't work superbly - it's kind of hard to do this in MATLAB afaik
% A possible better alternative would be to save the two arrays with
% georeferencing and then read them into ArcGIS or QGIS for comparison 

figure(3) 
imshow(RGB_im) 
hold on 

% plot NDWI areas as an overlay 
overlay1 = imagesc(NDWI_defined_SGLs);
ax1 = gca; 
colormap(ax1,hsv); 
newCData = overlay1.CData * 0.5; 
overlay1.AlphaData = newCData; 

%% Finish up 

cd(S2_main_dir); 
