% (c) Rebecca Dell, University of Cambridge, June 2020. 
% Adapted from the work of Williamson et al. (2018) (doi:10.5194/tc-2018-56). 

% Script to convert images from DN to TOA. 

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

% L-8 image bands were pre-downloaded from the USGS Earth Explorer interface
% Metadata files must also be named with same convention, ending in .txt. 

% Other scripts required = LS8_meta_parser, ls8_dn2refl_perpix. Created by
% Neil Arnold. See comments in these scripts for additional packages to
% download. 

clc; close all; clear all; 

%% Hard-code parameters and directories 

% change directory names according to the file locations on your computer
L8_main_dir = '/Users/rebeccadell/Desktop/l8_angles'; 
L8_in_folder = '/Users/rebeccadell/Desktop/l8_angles';
L8_out_folder = '/Users/rebeccadell/Desktop/l8_angles';

% shouldn't need to change the below unless you want the folder names to
% change... 
L8_out_bands_folder = [L8_out_folder '/bands'];

%% Load Landsat-8 data 

% make directories (if they don't already exist) 
mkdir(L8_out_folder); mkdir(L8_out_bands_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_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_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_b8; % checking number of b8 files 
check4 = 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 8 files does not match the number of band 2 files.'); 
elseif check4 ~= 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 % 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
    blue_band_fname = L8_b2_list(i).name; 
    green_band_fname = L8_b3_list(i).name; 
    red_band_fname = L8_b4_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
    info_geo_blue = geotiffinfo(blue_band_fname);
    [X_blue,R_blue] = geotiffread(blue_band_fname);
    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
    info_geo_pan = geotiffinfo(pan_band_fname);
    [X_pan,R_pan] = geotiffread(pan_band_fname);
%     geo_Tags_pan = info_geo_pan.GeoTIFFTags.GeoKeyDirectoryTag;
%     tiffTags = struct('TileLength',5255,'TileWidth',7653);
%     [B,S] = geotiffread(pan_band_fname); % reads in the image and stores the spatial data
    disp('Loading metadata file.') 
    meta = ls8_meta_parser(meta_fname); % Neil Arnold's function for L8
    
    % interpolate pan band to 30 m resolution (bilinearly following Pope
    % and others, 2016) 
    target_size = size(red_band); 
    pan_band = imresize(pan_band,target_size,'bilinear'); 
    
    % 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_perpix(imfilename,2); 
    red_band = ls8_dn2refl_perpix(imfilename,4); 
    green_band = ls8_dn2refl_perpix(imfilename,3); 
    pan_band = ls8_dn2refl_perpix(imfilename,8); 
        
     % 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'];
   
    % 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')
   geotiffwrite([new_blue_var_name '.tif'],blue_band,R_blue,'GeoKeyDirectoryTag', info_geo_blue.GeoTIFFTags.GeoKeyDirectoryTag)
    disp('Saving green band.') 
    save([new_green_var_name '.mat'],'green_band')
    geotiffwrite([new_green_var_name '.tif'],green_band,R_blue,'GeoKeyDirectoryTag', info_geo_blue.GeoTIFFTags.GeoKeyDirectoryTag)
    disp('Saving red band.') 
    save([new_red_var_name '.mat'],'red_band')
    geotiffwrite([new_red_var_name '.tif'],red_band,R_blue,'GeoKeyDirectoryTag', info_geo_blue.GeoTIFFTags.GeoKeyDirectoryTag)
 
    % N.B Here we save the panchromatic band back at a 15 m resolution. 
    disp('Saving panchromatic band.') 
    save([new_pan_var_name '.mat'],'pan_band', '-v7.3')
    geotiffwrite([new_pan_var_name '.tif'],pan_band,R_pan,'GeoKeyDirectoryTag', info_geo_pan.GeoTIFFTags.GeoKeyDirectoryTag)
 
    
    % 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 

cd(L8_main_dir); 

