% (c) Rebecca Dell, University of Cambridge, June 2020. 
% Adapted from the work of Williamson et al. (2018) (doi:10.5194/tc-2018-56). 
    % Who built on work of Pope et al. (2016)(doi:10.5194/tc-10-15-2016), Banwell et al. (2014)(doi:10.3189/2014AoG66A049), Arnold et al.
    % (2014)(doi:10.5194/tc-8-1149-2014), and Williamson et al. (2017) (doi:10.1016/j.rse.2017.04.032).

% Script to pre-process Landsat-8 data and produce a water body area and depth matrix for each image.   

% Landsat-8 (30 m or 15 m resolution) band specifications: 
% Band 2 = blue
% Band 3 = green 
% Band 4 = red
% Band 8 = panchromatic  

% L-8 image bands were pre-downloaded from the USGS Earth Explorer interface as DN values, converted to per-pixel TOA using the relevant scripts, 
% and then cropped and mosaicked to the ROI. 

% Bands have the same naming convention as the original downlaoded Landsat files (this must be adhered to 
% STRICTLY for this script to work). 

% For this study, cloud masks were made in GEE using 'Simple Cloud
% Score'(https://developers.google.com/earth-engine/landsat) and then
% downloaded for use in this script. 

% Other inputs include an ocean area mask and a study area mask
% Other scripts called within this script: S2_NDWI, L8_depth_mask

% Note whilst some variables refer to lakes (e.g. 'lake_depths_c'), at this
% stage they are just products showing 'all water bodies'. 

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 = 'C:\Users\rld46\Dropbox (Cambridge University)\testing_code'; 

if strcmp(year_spec,'2017') == 1   
    L8_in_folder = 'C:\Users\rld46\Dropbox (Cambridge University)\testing_code\input_L8';
elseif strcmp(year_spec,'2016') == 1 
    L8_in_folder = 'C:\Users\rld46\Dropbox (Cambridge University)\testing_code\input_L8';
end 
   
if strcmp(year_spec,'2017') == 1   
    L8_out_folder = 'C:\Users\rld46\Dropbox (Cambridge University)\testing_code\preprocessedL8';
elseif strcmp(year_spec,'2016') == 1 
    L8_out_folder = 'C:\Users\rld46\Dropbox (Cambridge University)\testing_code\preprocessedL8';
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'];


%% 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_b8_list = dir('*B8.tif'); % pan band list (must have .TIF extension) 
L8_cloud_mask_list =dir('*_CM.tif') % cloud mask list (.TIF 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_b8 = length(L8_b8_list); % number of band 8 files 
num_L8_files_cloud = length(L8_cloud_mask_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_b8; % checking number of b8 files 
check4 = num_L8_files_b2 - num_L8_files_cloud; % 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 8 files does not match the number of band 2 files.'); 
elseif check4 ~= 0 
    disp('Error! The number of cloud 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 check6 % no longer need these variables 

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 as it's nicer this way
    current_date = dates{i,1}; 
    year = current_date(1,20:21); %year_num = str2double(year); 
    month = current_date(1,22:23); %month_num = str2double(month); 
    day = current_date(1,24:25); %day_num = str2double(day);
    imfilename = current_date(1,1:40);
    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
    B2_fname = L8_b2_list(i).name; %Blue Band
    B3_fname = L8_b3_list(i).name; %Green Band
    B4_fname = L8_b4_list(i).name; %Red Band
    B8_fname = L8_b8_list(i).name; %Pan Band
    CM_fname = L8_cloud_mask_list(i).name; %Cloud Mask
    
    % read in the raw band files and assign variables 
    disp('Loading blue band data.') 
    B2 = imread(B2_fname); 
    % the two lines below also read in the georefferencing information
    info_geo = geotiffinfo(B2_fname);
    [X,R] = geotiffread(B2_fname); 
    B2_forRGB = B2; 
    B2 = double(B2); % convert to double for processing
    disp('Loading green band data.') 
    B3 = imread(B3_fname); 
    B3_forRGB = B3; 
    B3 = double(B3); % convert to double for processing
    disp('Loading red band data.') 
    B4 = imread(B4_fname); 
    B4_forRGB = B4; 
    B4 = double(B4); % convert to double for processing
    disp('Loading panchromatic band data.') 
    B8 = imread(B8_fname); 
    B8 = double(B8); % convert to double for processing
    disp('Loading cloud mask data.') 
    CM = imread(CM_fname); 
    CM = double(CM); % convert to double for processing
    
    % interpolate pan band to 30 m resolution (bilinearly following Pope
    % and others, 2016) 
    target_size = size(B4); 
    B8 = imresize(B8,target_size,'bilinear'); 
      
    % make a mask to show the locations that have no data for them 
    no_data_mask = zeros(size(B4)); 
    no_data_mask(B4 <= 0) = 1; % all negative values in L8 are no data
    
    % mask no-data values (negatives in LS8) 
    B8(B8 <= 0) = NaN; 
    B2(B2 <= 0) = NaN; 
    B4(B4 <= 0) = NaN;
    B3(B3 <= 0) = NaN;
    
    % mask cloudy pixels 
    B8(CM ==1) = NaN; 
    B2(CM ==1) = NaN; 
    B4(CM ==1) = NaN;
    B3(CM ==1) = NaN;
   
    %% Create RGB Image
    
    % show RGB image and stretch to the limits of the data 
    RGB_out = cat(3,B4_forRGB,B3_forRGB,B2_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')
    
%% Set Constant NDWI testing values 

LS8_pix_siz = 0.03; LS8_pix_siz_m = 30; % specifying L8 resolution 
cd(L8_in_folder);
load('test21_S2.mat'); %ocean test area mask (1's on ocean)
load('ice_ocean_mask_S2.mat')%study area mask (1's on ice)

%Resample masks from S2 resolution to L8 resolution 
target_size = size(B4);
test21_S2 = imresize(test21_S2,target_size,'nearest'); 
ice_ocean_mask_S2 = imresize(ice_ocean_mask_S2,target_size,'nearest'); 
    
g_red = 0.7507; % value from Pope and others (2016) 
g_pan = 0.3817; % value from Pope and others (2016) 

test21_S2 = double(test21_S2);
tminf_pan = B8; %makes a new matrix identical to the pan band for getting rinf values
tminf_red = B4; %makes a new matrix identical to the red band for getting rinf values
tminf_pan(test21_S2 == 0) = NaN; % test21 is a mask of the ocean, which ensures the inf values calculated are taken from optically deep water and not water on the nunatak.
tminf_red(test21_S2 == 0) = NaN;
% get darkest pixel values for the images for calibrating the
% depth-reflectance curve (see Williamson and others, 2017, for
% details, based on Sneed and Hamilton, 2007). 
rinf_pan = min(min(tminf_pan)); 
rinf_red = min(min(tminf_red)); 
rinf_pan_str = num2str(rinf_pan,'%.3f'); 
rinf_red_str = num2str(rinf_red,'%.3f'); 
disp(['Darkest red band pixel value is ' rinf_red_str '.']); 
disp(['Darkest panchromatic band value is ' rinf_pan_str '.']);
    
      if rinf_red > 0.1 
        disp('Warning: darkest red band value is high - manually check.');
        rinf_red = 0.0    
      elseif rinf_pan > 0.1 
        disp('Warning: darkest panchromatic band value is high - manually check.'); 
        rinf_pan = 0.0 
      end

%% NDWI Area Calculations
% classify SGLs for this image with the NDWI technique  

n = 1

% Set the NDWI threshold
x = 0.25 % repeat for each file in list 
    
%   disp(['Now analysing threshold' x;
    NDWI_str = num2str(x); 
    disp(['Classifying wet pixels based on the NDWI technique (threshold = ' NDWI_str ').'])  
    wet_pixels = S2_NDWI(B2,B4,x); % NDWI function (another script called in here for this).
    
% Make a water body matrix
    cc = bwconncomp(wet_pixels,8); %Identifies 'objects' with '8' connectivity
    stats = regionprops(cc,'Area'); % Find the area of each object
    water_bodies_find = find([stats.Area] > 2); % From cc, finds objects with an area over two pixels
    water_bodies_stats = stats(water_bodies_find); % Creates a new table/ variable, which just has the objects identified in the line above
    combined_matrix = ismember(labelmatrix(cc),water_bodies_find); % Creates a new matrix from identified objects
    combined_matrix= double(combined_matrix);%converts the variable from a logical to a double
    combined_matrix_filled = imfill(combined_matrix,'holes'); % Fill in the holes of each water body
    
    % Add the non-filled and filled matricies to create a mask of 'filled
    % pixels' dry = 0, filled =1, always wet =2. Get rid of always wet so we
    % are just looking at the pixels that have been filled. This will be needed
    % for the depth calculations.
    filled_combined = combined_matrix + combined_matrix_filled;
    filled_combined(filled_combined == 2) = 0;
    
%% Depth Calculations
 
    combined_matrix(isnan(combined_matrix)) = 0; % converts any NaNs in the water body mask to 0
    combined_matrix(ice_ocean_mask_S2 == 0) = 0;
    combined_matrix_filled(isnan(combined_matrix_filled)) = 0; % converts any NaNs in the water body mask to 0
    combined_matrix_filled(ice_ocean_mask_S2 == 0) = 0;
    filled_combined(isnan(filled_combined)) = 0; % converts any NaNs in the water body mask to 0
    filled_combined(ice_ocean_mask_S2 == 0) = 0;
    
    % This is the depth function
   [lake_depths_c] = L8_depth_mask(rinf_red,rinf_pan,B8,B4,combined_matrix,filled_combined);
    
   lake_depths_c(lake_depths_c == 0) = NaN;

%% Make 10m masks for comparison with Sentinel

    % resample combined mask to 10 m for water body tracking process 
    size_30m = size(combined_matrix_filled); % get the 30 m size 
    target_size = size_30m .* 3; % multiply by 3 to get to 10 m resolution
    combined_matrix_filled_10m = imresize(combined_matrix_filled,target_size,'nearest'); %make 10m water body mask
    lake_depths_c_10m = imresize(lake_depths_c,target_size,'bilinear'); % make 10m depth matrix
    
% Save Variables

    % set new variable names with the short date string 
    new_blue_var_name = [imfilename '_B2']; 
    new_green_var_name = [imfilename '_B3']; 
    new_red_var_name = [imfilename '_B4']; 
    new_pan_var_name = [imfilename '_B8'];
    new_RGB_var_name = [imfilename '_RGB'];
    new_cmask_var_name = [imfilename 'cloud_mask'];
    new_no_data_var_name = [imfilename '_no_data_mask'];
    new_combined_matrix_name = [imfilename '_combined_matrix'];
    new_combined_matrix_10m_name = [imfilename '_combined_matrix_10m'];
    new_depth_mask_name = [imfilename '_depths'];
    new_depth_mask_name_10m = [imfilename '_depths_10m'];
    
    %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'],'B2')
    disp('Saving green band.') 
    save([new_green_var_name '.mat'],'B3')
    disp('Saving red band.') 
    save([new_red_var_name '.mat'],'B4')
    disp('Saving SWIR band.') 
    save([new_pan_var_name '.mat'],'B8')
    cd(L8_out_NDWI_folder); 
    disp('Saving combined_matrix.')
    save([new_combined_matrix_name '.mat'],'combined_matrix_filled')
    disp('Saving lake_depths_c.')
    save([new_depth_mask_name '.mat'],'lake_depths_c')
    
    
    cd(L8_out_NDWI_10m_folder);
    disp('Saving combined_matrix_filled_10m.')
    save([new_combined_matrix_name '.mat'],'combined_matrix_filled_10m')
    disp('Saving lake_depths_c_10m.')
    save([new_depth_mask_name_10m '.mat'],'lake_depths_c_10m')

    
% Saving the water mask and the depth matrix as geotiffs.   
  geotiffwrite([new_depth_mask_name '.tif'],lake_depths_c,R,'GeoKeyDirectoryTag', info_geo.GeoTIFFTags.GeoKeyDirectoryTag)
  cd(L8_out_NDWI_folder); 
  geotiffwrite([new_combined_matrix_name '.tif'],combined_matrix_filled,R,'GeoKeyDirectoryTag', info_geo.GeoTIFFTags.GeoKeyDirectoryTag)  

   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'],'CM')
   cd(L8_out_no_data_folder); 
   disp('Saving No Data mask.') 
   save([new_no_data_var_name '.mat'],'no_data_mask'); 

    % 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(' ') 
    
    
    
end 
 
