%% (c) Rebecca Dell, University of Cambridge, June 2020
% Script to split the full water body depth array's into two arrays (representing linear and circular water bodies) based on a solidity threshold

% Note that circular water bodies are called lakes throughout the script,
% and linear water bodies are called streams. This terminology is
% incorrect, but was simply used at the time of writing the code.  

%% Hard-code parameters and directories 
% Adjust the file paths to either work on landsat or sentinel data-sets as
% appropriate
root_dir = 'C:\Users\rld46\Dropbox (Cambridge University)'; 
depth_dir = 'C:\Users\rld46\Dropbox (Cambridge University)\testing_code\preprocessedS2\NDWI_masks'; % or preprocessedL8\NDWI_10m_masks
out_dir = 'C:\Users\rld46\Dropbox (Cambridge University)\testing_code\CategorisedS2'; %or CategorisedL8

% find depths 
cd(depth_dir); 
depth_list = dir('*_depths.mat'); % OR ('*_depths_10m.mat');
 
% % get the number of images from the length of the depth list 
img_num = length(depth_list);

% set up a blank array with same number of rows as img_num and 3 columns 
date_nums = zeros(img_num,3); 
dates = {depth_list(:).name}; % get the names of the images
dates = dates'; % make vertical 

% for each of the dates in list, do this... 
for i = 1:img_num 
    
    current_date = dates{i,1}; 
    
        if current_date(1,1) == 'L'
        image_type = 'L8'
        year = current_date(1,20:21);
        month = current_date(1,22:23); 
        day = current_date(1,24:25);
        imfilename = current_date(1,20:25);
    else
        image_type = 'S2'
        year = current_date(1,1:2);
        month = current_date(1,3:4);
        day = current_date(1,5:6);
        imfilename = current_date(1,1:6);
        end
    
    curr_date_str_shrt = [year month day]; 
    curr_date_num = str2double(curr_date_str_shrt); % convert to double 
    date_nums(i,2) = curr_date_num; % put the date into 2nd column 
    date_nums(i,1) = i; % put the number into 1st column
    
    curr_depth_array_name = depth_list(i).name;
    load(curr_depth_array_name); 

    curr_depth_array = lake_depths_c; %(SENTINEL) for Landsat add _10m
    curr_area_array= zeros(size(curr_depth_array));
    curr_area_array(curr_depth_array > 0) = 1;

    % Create a binary image with an area over X and an Solidity Index less than Y    
    cc = bwconncomp(curr_area_array,8); %Identifies 'objects' with '8' connectivity
    stats = regionprops(cc,'Solidity');
    lakes_find = find([stats.Solidity] >= 0.45); % From cc,finds objects with a solidity >= 0.45
    lakes_stats = stats(lakes_find); % Creates a new table/ variable, which just has the objects identified in the previous stage
    curr_lake_area_array = ismember(labelmatrix(cc),lakes_find); % Creates a new matrix from identified objects
    curr_lake_area_array= double(curr_lake_area_array); %converts the variable from a logical to a double

    % Create a binary image with an Eccentricity Index less
    % than Y (Same steps as above, just for streams)
    streams_find = find([stats.Solidity] < 0.45);  
    streams_stats = stats(streams_find);
    curr_stream_area_array = ismember(labelmatrix(cc),streams_find); 
    curr_stream_area_array = double(curr_stream_area_array);
    curr_stream_area_array(curr_stream_area_array > 0) = 2; % changes all 'stream pixels' to a value of 2, so we can later combine the lakes and streams matricies. 

    curr_lake_depth_array = curr_depth_array;
    curr_lake_depth_array(curr_lake_area_array == 0) =0;
    curr_stream_depth_array = curr_depth_array;
    curr_stream_depth_array(curr_stream_area_array == 0) =0;

    curr_combined_area_array = curr_lake_area_array + curr_stream_area_array;
    curr_combined_depth_array = curr_lake_depth_array + curr_stream_depth_array;

    new_lake_depth_mask_name = [image_type imfilename '_depths_l']; 
    new_stream_depth_mask_name = [image_type imfilename '_depths_s'];
    new_lake_area_mask_name = [image_type imfilename '_area_l'];
    new_stream_area_mask_name = [image_type imfilename '_area_s'];
    new_combined_depth_mask_name = [image_type imfilename '_depths_c'];
    new_combined_area_mask_name = [image_type imfilename '_area_c'];

    %save the bands to the output folders and remove from workspace 
    cd(out_dir);
    disp('Saving depths') 
    save([new_lake_depth_mask_name '.mat'],'curr_lake_depth_array')
    save([new_stream_depth_mask_name '.mat'],'curr_stream_depth_array')
    save([new_combined_depth_mask_name '.mat'],'curr_combined_depth_array')
    disp('Saving areas') 
    save([new_lake_area_mask_name '.mat'],'curr_lake_area_array')
    save([new_stream_area_mask_name '.mat'],'curr_stream_area_array')
    save([new_combined_area_mask_name '.mat'],'curr_combined_area_array')
    cd(depth_dir);
        
end