
% Script to pre-process Landsat-8 data. 

% Landsat-8 (30 m or 15 m resolution) band specifications: 
% Band 2 = blue
% Band 3 = green 
% Band 4 = red
% Band 6 = SWIR (for clouds) 
% Band 8 = panchromatic (for SGL depths) 

% L-8 image bands were pre-downloaded from the USGS Earth Explorer interface, 
% then cropped to the same extent of the S-2 images in ArcGIS. 
% Bands have the following naming convention (this must be adhered to 
% STRICTLY for this script to work): 
% YEAR-MO-DY-BN.TIF, 
% where YEAR = year (e.g. 2016), MO = month (e.g. 08), DY = day (e.g. 31)
% and N = band number (e.g. 2). Example file name: 2016-05-15-B3.TIF 
% Metadata files must also be named with same convention, ending in .txt. 

% Andrew Williamson, University of Cambridge, July 2017 

clc; close all; clear all; 

%% Hard-code parameters and directories 

year_spec = '2017'; 

% change directory names according to the file locations on your computer
L8_main_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/'; 

if strcmp(year_spec,'2017') == 1   
    L8_in_folder = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/raw_LS8_images';
elseif strcmp(year_spec,'2016') == 1 
    L8_in_folder = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/raw_LS8_images';
end 
   
if strcmp(year_spec,'2017') == 1   
    L8_out_folder = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_LS8';
elseif strcmp(year_spec,'2016') == 1 
    L8_out_folder = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_LS8';
end 

% shouldn't need to change the below unless you want the folder names to
% change... 
L8_out_bands_folder = [L8_out_folder '/bands'];
L8_out_RGBs_folder = [L8_out_folder '/RGBs'];
L8_out_RGB_images_folder = [L8_out_folder '/RGB_images'];
L8_out_NDWI_folder = [L8_out_folder '/NDWI_masks'];
L8_out_cmask_folder = [L8_out_folder '/cloud_masks'];
L8_out_NDWI_10m_folder = [L8_out_folder '/NDWI_10m_masks'];
L8_out_no_data_folder = [L8_out_folder '/no_data_masks'];

NDWI_threshold = 0.25; % set the NDWI threshold for classifying SGLs 
% recommend testing a variety of thresholds before incorporating here 
save_to_wkspace = 'N'; % save variables to workspace? 'Y' or 'N'
% if specifying 'Y' and there are lots of files, then MATLAB 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 

% note that cloud masks here are saved at 30 m resolution, so I resample
% them later in the lake tracking to 10 m for using them within the lake
% tracking (using a nearest-neighbour resampling regime) 

%% Load Landsat-8 data 

% make directories (if they don't already exist) 
mkdir(L8_out_folder); mkdir(L8_out_bands_folder); mkdir(L8_out_RGBs_folder);
mkdir(L8_out_NDWI_folder); mkdir(L8_out_RGB_images_folder); 
mkdir(L8_out_cmask_folder); mkdir(L8_out_NDWI_10m_folder); 
mkdir(L8_out_no_data_folder); 
addpath(L8_main_dir); 
cd(L8_in_folder); % change to input folder 
L8_b2_list = dir('*B2.TIF'); % blue band list (must have .TIF extension) 
L8_b3_list = dir('*B3.TIF'); % green band list (must have .TIF extension) 
L8_b4_list = dir('*B4.TIF'); % red band list (must have .TIF extension) 
L8_b6_list = dir('*B6.TIF'); % swir band list (must have .TIF extension) 
L8_b8_list = dir('*B8.TIF'); % pan band list (must have .TIF extension) 
L8_meta_list = dir('*MTL.txt'); % metadata list (.txt extension) 
num_L8_files_b2 = length(L8_b2_list); % number of band 2 files 
num_L8_files_b3 = length(L8_b3_list); % number of band 3 files 
num_L8_files_b4 = length(L8_b4_list); % number of band 4 files 
num_L8_files_b6 = length(L8_b6_list); % number of band 6 files 
num_L8_files_b8 = length(L8_b8_list); % number of band 8 files 
num_L8_files_meta = length(L8_meta_list); % number of meta data files 
string = num2str(num_L8_files_b2); % convert to string for command line 

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

check1 = num_L8_files_b2 - num_L8_files_b3; % checking number of b3 files 
check2 = num_L8_files_b2 - num_L8_files_b4; % checking number of b4 files 
check3 = num_L8_files_b2 - num_L8_files_b6; % checking number of b6 files 
check4 = num_L8_files_b2 - num_L8_files_b8; % checking number of b8 files 
check5 = num_L8_files_b2 - num_L8_files_meta; % check number metadata 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.'); 
elseif check3 ~= 0 
    disp('Error! The number of band 6 files does not match the number of band 2 files.'); 
elseif check4 ~= 0 
    disp('Error! The number of band 8 files does not match the number of band 2 files.'); 
elseif check5 ~= 0 
    disp('Error! The number of metadata 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 check3 check4 check5 % no longer need these variables

% no-data and cloud info arrays set-up 
no_data_info = zeros(1,4); % zeros array with 1 row and 4 columns 
cloud_info = zeros(1,4); % zeros array with 1 row and 4 columns 
    
for i = 1:num_L8_files_b2 % repeat for each file in list 
     
    % first get the date information from the input file
    % only need to do it for band 2 since it should match for all
    % anyway, assuming they are in the file lists in the same order 
    dates = {L8_b2_list(:).name}; 
    dates = dates'; % transpose matrix 
    current_date = dates{i,1}; 
    year = current_date(1,3:4); %year_num = str2double(year); 
    month = current_date(1,6:7); %month_num = str2double(month); 
    day = current_date(1,9:10); %day_num = str2double(day); 
    curr_date_str = [day '/' month '/' year]; % format date string nicely 
    curr_imnum_str = num2str(i); 
    tot_im_str = num2str(num_L8_files_b2); 
    
    disp(['Now analysing image from ' curr_date_str ', image number '...
        curr_imnum_str ' of ' tot_im_str '.']); 
    
    % new variable for short name string
    short_dt_str = ['20' year month day]; 
    
    % find out which variable to read in
    blue_band_fname = L8_b2_list(i).name; 
    green_band_fname = L8_b3_list(i).name; 
    red_band_fname = L8_b4_list(i).name; 
    SWIR_band_fname = L8_b6_list(i).name; 
    pan_band_fname = L8_b8_list(i).name;
    meta_fname = L8_meta_list(i).name; 
    
    % read in the raw band files and assign variables 
    cd(L8_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
    disp('Loading panchromatic band data.') 
    pan_band = imread(pan_band_fname); 
    pan_band = double(pan_band); % convert to double for processing
    disp('Loading SWIR band data.') 
    SWIR_band = imread(SWIR_band_fname); 
    SWIR_band = double(SWIR_band); % convert to double for processing
    disp('Loading metadata file.') 
    meta = ls8_meta_parser(meta_fname); % Neil Arnold's function for L8
    
    % convert digital numbers to radiance for the bands 
    % see the function script for details; again, Neil Arnold's function. 
    % the numbers specify the original band of L8
    blue_band = ls8_dn2refl(blue_band,meta,2); 
    red_band = ls8_dn2refl(red_band,meta,4); 
    green_band = ls8_dn2refl(green_band,meta,3); 
    SWIR_band = ls8_dn2refl(SWIR_band,meta,6); 
    pan_band = ls8_dn2refl(pan_band,meta,8); 
    
    % create cloud mask where band-6 reflectance is > 0.1 (following
    % Williamson et al. (2017) but with a more stringent threshold of 0.1,
    % not 0.15 as with MODIS) 
    cloud_mask = zeros(size(SWIR_band)); % 0s array  with size of SWIR band
    cloud_mask(SWIR_band > 0.1) = 1; % set as 1s if > 0.1 
    se = strel('disk',2); 
    cloud_mask = imdilate(cloud_mask,se); 
    
    % calculate the percentage of cloud cover relative to whole image 
    total_cloud = sum(sum(cloud_mask)); 
    row_num = size(cloud_mask,1); 
    col_num = size(cloud_mask,2); 
    image_size = row_num * col_num; 
    pct_clouds = 100/image_size * total_cloud;
    cloud_info(i,1) = i; % image number 
    cloud_info(i,2) = image_size; % total number of pixels 
    cloud_info(i,3) = total_cloud; % total number of no-data pixels 
    cloud_info(i,4) = pct_clouds; % percentage no-data pixels
    
    % make a mask to show the locations that have no data for them 
    no_data_mask = zeros(size(red_band)); % 0s array  with size of red band
    no_data_mask(red_band < 0) = 1; % all negative values in L8 are no data
    
    % calculate the percentage of no-data values relative to whole image 
    total_no_data = sum(sum(no_data_mask)); 
    row_num = size(no_data_mask,1); 
    col_num = size(no_data_mask,2); 
    image_size = row_num * col_num; 
    pct_no_data = 100/image_size * total_no_data; 
    pct_no_data = 100 - pct_no_data; 
    no_data_info(i,1) = i; % image number 
    no_data_info(i,2) = image_size; % total number of pixels 
    no_data_info(i,3) = total_no_data; % total number of no-data pixels 
    no_data_info(i,4) = pct_no_data; % percentage no-data pixels 
    
    %%
    
    % 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(L8_out_RGB_images_folder);
    disp('Saving RGB image.') 
    out_im_name = ['L8_' 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  
    
    % mark-cloud masked lakes as NaNs 
    NDWI_defined_SGLs(cloud_mask == 1) = NaN; 
    
    % mark no-data values as NaNs 
    NDWI_defined_SGLs(no_data_mask == 1) = NaN; 
    
    % resample NDWI mask to 10 m for lake tracking process 
    size_30m = size(cloud_mask); % get the 30 m size from cloud_mask size
    target_size = size_30m .* 3; % multiply by 3 to get to 10 m resolution
    L8_lmask_10m = imresize(NDWI_defined_SGLs,target_size,'nearest'); 
    % make new lake mask at 10m using these parameters and 'nearest
    % neighbour' to avoid changing any of the original values (which
    % another interpolation regime, like bilinear, would involve) 
    
    % filter small lakes (< 5 pixels) and linear lakes with function 
    % the 2 here refers to the linear features too (< 2 pixels wide) 
    NDWI_defined_SGLs = SGL_mask_filter(NDWI_defined_SGLs,5,2);
    
    % set new variable names with the short date string 
    new_blue_var_name = ['L8_' short_dt_str '_blue']; 
    new_green_var_name = ['L8_' short_dt_str '_green']; 
    new_red_var_name = ['L8_' short_dt_str '_red'];
    new_SWIR_var_name = ['L8_' short_dt_str '_SWIR'];
    new_pan_var_name = ['L8_' short_dt_str '_pan'];
    new_RGB_var_name = ['L8_' short_dt_str '_RGB'];
    new_NDWI_var_name = ['L8_' short_dt_str '_NDWI_SGL_mask'];
    new_NDWI_10m_var_name = ['L8_' short_dt_str '_NDWI_10m_SGL_mask'];
    new_cmask_var_name = ['L8_' short_dt_str '_cloud_mask'];
    new_no_data_var_name = ['L8_' short_dt_str '_no_data_mask'];
    
    % save the bands to the output folders and remove from workspace 
    cd(L8_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')
    disp('Saving SWIR band.') 
    save([new_SWIR_var_name '.mat'],'SWIR_band')
    disp('Saving panchromatic band.') 
    save([new_pan_var_name '.mat'],'pan_band')
    cd(L8_out_NDWI_folder); 
    disp('Saving NDWI mask.') 
    save([new_NDWI_var_name '.mat'],'NDWI_defined_SGLs')
    cd(L8_out_RGBs_folder); 
    disp('Saving RGB array.') 
    save([new_RGB_var_name '.mat'],'new_RGB_image')
    cd(L8_out_cmask_folder); 
    disp('Saving cloud mask.') 
    save([new_cmask_var_name '.mat'],'cloud_mask')
    disp('Saving NDWI mask at 10 m resolution.')
    cd(L8_out_NDWI_10m_folder); 
    save([new_NDWI_10m_var_name '.mat'],'L8_lmask_10m')
    disp('Saving no-data mask at 10 m resolution.') 
    no_data_mask_10m = imresize(no_data_mask,target_size,'nearest'); 
    cd(L8_out_no_data_folder); 
    save([new_no_data_var_name '.mat'],'no_data_mask_10m'); 
    
    % 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_pan_var_name,pan_band);
        assignin('base',new_SWIR_var_name,SWIR_band);
        assignin('base',new_RGB_var_name,new_RGB_image);
        assignin('base',new_NDWI_var_name,NDWI_defined_SGLs);
        assignin('base',new_cmask_var_name,cloud_mask);
        assignin('base',new_NDWI_10m_var_name,L8_lmask_10m);
        assignin('base',new_no_data_var_name,no_data_mask_10m);
    end
    
    % get percentage progress update for command line and post to it 
    pct_complete = (100/num_L8_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 pan_band ...
        SWIR_band cloud_mask size_30m L8_lmask_10m target_size ...
        new_NDWI_10m_var_name no_data_mask meta_fname meta ...
        new_cmask_var_name new_pan_var_name new_SWIR_var_name ...
        out_im_name pan_band_fname SWIR_band_fname ...
        total_no_data row_num col_num image_size pct_no_data ...
        total_cloud pct_clouds se curr_imnum_str no_data_mask_10m ...
        new no_data_var_name 
    % keeping workspace tidy (and small enough to allow processing!)
    end 
    
end 

cd(L8_main_dir); 
