
% Script to ascertain rapidly and slowing draining SGLs from a dataset
% based upon defined parameters. This is setup for Sentinel-2 and Landsat 8
% lake tracking but could be adjusted to other sets of imagery. 

% While "slowly" draining lakes are identified here, in reality every lake
% in the dataset that is not identified as "rapidly" draining is identified
% as draining slowly because they all lose > 20% of their volume over any
% time period. In Miles et al. (2017, Front. Earth Sci.), they 
% further separated these lakes into slowly draining and freezing by using
% a manual method. Without having access to e.g. climate data, this cannot
% be done. This is a possible future improvement to this work. 

% Can save outputs as Excel files if wished, but I've commented this out
% because I don't like Excel... 

% Andrew Williamson, University of Cambridge, January 2018 

%% Specify input parameters for running 

imagery = 'Sentinel'; % specify here based on which imagery type you're 
% using based on which one the lake tracking was done for, either
% 'Landsat','Sentinel' or 'both' 

drainage_type = 'volume'; % will drainage events be identified based on 
% area or volume changes?  Note that even if you choose 'volume' here,
% lakes are still separated into 'large' and 'small' based on their area
% arrays (as that's the threshold that affects MODIS identification, not
% the volume) 

rpd_loss_reqd = 80; % whats percent of volume or area is required to define 
% a rapid event? Note that this is 80 % of the lake's ** maximum seasonal
% value ** and not an individual daily value

slw_loss_reqd = 20; % same idea for slow lake drainage, specifying 20 % 
% here, which is of the ** maximum seasonal value ** 

rapid_time_d = 4; % period (in days) over which rapid lake drainage is 
% identified  

allow_refill = 'Y'; % for rapid events, do you want to allow refilling of 
% the lake basin following the classification of a drainage event and for
% the drainage event to still be classified as having happened? choose
% either 'Y' if you want to allow this or 'N' if you do not, instead 
% wanting the basin to remain completely empty

refill_amt_to_allow = 20; % if 'Y' above, how much refilling to allow
% (in percent) compared with the starting value? rationale behind this 
% 20% figure is to ensure that even with refilling to this value, it would
% still have been classified as a lake-drainage event since it would still
% have lost >= 80% of its starting area. N.B. this is for rapid events! 

allow_slw_refill = 'Y'; % asking same question regarding refilling as above 
rfill_slw_amt_allow = 80; % and by how much if 'Y'? 

slow_time_d = 1; % slow drainage of this sort can occur over any time period

tracking_area_big = 1250; % area of 'big' lakes that are normally possible
% to be tracked with MODIS imagery, which are outputted differently 

pix_size = 0.01; % Sentinel-2 pixel size is 0.01 km (10 m) 

post_drain_cloudy_days = 1; % number of post-drainage days to take into
% account when examining refilling of lakes and ensuring that a basin
% remains empty following drainage. Deliberately using a less stringent 
% threshold here because the time interval between images for this dataset
% is relatively high compared with e.g. MODIS. 

tidy_wkspace = 'N'; % tidy up workspace to remove redundant variables? 'Y'
% or 'N' 

%% Specify and make directories 

disp(['Doing lake drainage analysis for ' imagery '.']) 
inputs = ['/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_lake_tracking_2017/' imagery '/']; 
% points to the correct directory for the inputs depending upon which has
% been used 
outputs = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_lake_drainage_2017'; 
mkdir(outputs); 
background_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/'; % main directory 
final_out_dir = [outputs '/' imagery '/' drainage_type '/']; 
mkdir(final_out_dir); 
all_outputs = [final_out_dir '/all_SGLs']; % folder for all lakes together
big_outputs = [final_out_dir '/big_SGLs']; % folder for just big lakes
small_outputs = [final_out_dir '/small_SGLs']; % folder for just small lakes 
mkdir(all_outputs); mkdir(big_outputs); mkdir(small_outputs); 

%% Load input files 

cd(inputs) 
% load the correct tracking array based on whether identifying lake
% drainage based on area or volume changes to lakes 
if strcmp(drainage_type,'area') == 1 
    
    % get correct input file name 
    fname = [inputs 'lake_area_tracking.mat']; 
    load(fname); 
    sampleSGLs = area_info_file; % specify variable accordingly 
    final_row_num = size(area_info_file,1); % identify the final row as it
    % contains the image dates, from which we obviously don't want to
    % identify lake drainage events... 
    dates = sampleSGLs(final_row_num,:); % get the dates from this row 
    sampleSGLs(final_row_num,:) = []; % delete this row of the array 
    
elseif strcmp(drainage_type,'volume') == 1 
    
    % get correct input file name 
    fname = [inputs 'lake_volume_tracking.mat'];   
    load(fname); 
    sampleSGLs = volume_info_file; % specify variable accordingly 
    final_row_num = size(volume_info_file,1); % identify final row as it
    % contains the image dates, from which we obviously don't want to
    % identify lake drainage events... 
    dates = sampleSGLs(final_row_num,:); % get the dates from this row 
    sampleSGLs(final_row_num,:) = []; % delete this row of the array 
    
end 

lake_number = size(sampleSGLs,1); % number of lakes is first dimension
image_number = size(sampleSGLs,2); % number of images is second dimension

rapid_time_h = rapid_time_d * 24; % gets drainage time in days to hours 
slow_time_h = slow_time_d * 24; % get drainage time in days to hours 

%% Convert dates (format in which all files were saved) to days of year 

% note that format of date was of the type e.g. 20160529, i.e. YYYYMMDD

all_JD_nums = zeros(size(dates,2),1); % new array with right sizes 

% repeat for-loop for the size of the dates array 
for i = 1:size(dates,2) 
    
    curr_date = dates(1,i); % take from i'th column 
    curr_date = num2str(curr_date); % convert to string 
    day = curr_date(7:8); % extract correct parts of string for day
    month = curr_date(5:6); % same for month 
    year = curr_date(1:4); % and year 
    d = [month '/' day '/' year]; % put in slashes (MATLAB date format) 
    v = datevec(d); % make this string into a date vector
    v0 = v; % copy this vector 
    v0(:,2:3) = 1; % take away the year from the data vector 
    DOY = datenum(v) - datenum(v0) + 1; % get DOY using this (accounts
    % for leap year if relevant) 
    all_JD_nums(i,1) = DOY; % put this info into the all_JD_nums array 
    
end 

%% Start main script functions 

current_events_row = 0; % this needs setting at 0 outside of the for-loop
% as it increments by 1 within it once an event has been identified 

drainage_events = zeros(1,10); 
% pre-allocate a base array of this size to speed up script 

for i = 1:lake_number % do this for each of the lakes within the max extent 
    
    current_SGL_values = sampleSGLs(i,:); 
    % get the current lake's values by extracting every column of current
    % row
    
    max_seasonal_SGL_area = nanmax(current_SGL_values); 
    % find out the maximum area/volume of the SGL over the season (using 
    % the "max" function may result in different output to "nanmax" here)
    % note that even though it says "area" here, if you choose "volume" in
    % the options at the start of the script, it will take into account
    % volume and not area 
    
    critical_loss = max_seasonal_SGL_area * (rpd_loss_reqd / 100); 
    % define the critical loss required to define a rapid SGL event based
    % upon the critical threshold inputted to the script and multiplying
    % the maximum seasonal area by this value 
    
    for p = 1:(image_number - 1) % repeat this for-loop for every image but
        % don't do it for the last one as obviously there are no data
        % beyond it at that point to compare it with...
        
        current_val = current_SGL_values(1,p); 
        % get the current SGL value by defining it using the "p" variable 
        
        next_values = current_SGL_values(1,(p + 1):end); 
        % extract everything from the current value + 1 to the end of the
        % array, meaning if the value is towards the start of season, there
        % are more data points against which this will be compared and so
        % the size of this array will vary over time as the for-loop
        % executes 
        
        number_nexts = size(next_values,2); 
        % find out how many "next_values" there are to compare the current
        % value against 
        
        for h = 1:number_nexts % do the following for each of the total
            % number of next values 
            
            next_val = next_values(1,h); 
            % get the current value being examined 
            
            if next_val < (current_val - critical_loss) 
                % compare the current value against the next value and take
                % the difference between them; if the next_value is less
                % than the current_value minus the critical loss required
                % for a drainage event, then do the following 
                
                % write the following out to the command line - note that
                % it's only a preliminary event because these need tidying
                % up to avoid duplicates and so on (as done below) 
%                 disp(['Found preliminary lake drainage event for ' ...
%                     'lake number ' num2str(i) ' beginning on image number '...
%                     num2str(p) '.']); 
                
                % now put all of the relevant information relating to this
                % event into the drainage_events array: 
                current_events_row = current_events_row + 1; 
                % increment this variable by 1 to make sure current data
                % are not being overwritten obvs 
                
               	drainage_events(current_events_row,1) = i; 
                % put the SGL number into the first column 
                start_JDay = all_JD_nums(p,1); % get the JD number for 
                % this day (although this is not used here - the JD numbers
                % are added to this array later for various reasons and 
                % replace the current day numbers) 
                
                drainage_events(current_events_row,2) = p;
                % put the current start day into column 2 
                
                % get the same information for the end day 
                end_day = h + p; 
                end_JDay = all_JD_nums(end_day,1); 
                drainage_events(current_events_row,3) = end_day; 
                % put this information into row 3 
                
                drainage_events(current_events_row,4) = current_val; 
                % find out the current value being examined on the start
                % day and put this into column 4 
                
                drainage_events(current_events_row,5) = next_val; 
                % same for end area, put into column 5 
                
                change = current_val - next_val; 
                drainage_events(current_events_row,6) = change; 
                % get the change in area and put this into column 6,
                % meaning that column 4 minus column 5 should equal column
                % 6 always
                
                percent_loss = 100 - (100 / current_val * next_val); 
                drainage_events(current_events_row,7) = percent_loss; 
                % also get the percentage loss to ensure that events are
                % not being identified badly for whatever reason - put this
                % information into column 7 
                
                change_period = (end_JDay - start_JDay) * 24; 
                drainage_events(current_events_row,8) = change_period;
                % work out period of time over which this event occurred by
                % taking end from start day and multiplying by 24 hours... 
                
                mean_drain_date = (start_JDay + end_JDay) / 2; 
                drainage_events(current_events_row,9) = mean_drain_date;
                % mean drainage date (I add JD number here at this stage as
                % this information is not used later) is calculated and put
                % into column 9 
                
                drain_date_error = (end_JDay - start_JDay) / 2; 
                drainage_events(current_events_row,10) = drain_date_error;
                % error in drainage date is in column 10 
                
            end % end of if next_value < critical loss test          
        end % end of for each of the "next_values" 
    end % end of do it for each image in the dataset 
end % end of for each lake number in the dataset 
    
number_events = size(drainage_events,1);
% work out how many events there are in the array at the moment by taking
% the size along the first (i.e. rows) dimension 

%% Filter events where basin refills above critical threshold post-drainage

% repeat the following loop for each of the number_events identified in the
% array in total 
for i = 1:number_events 
    
    contem_day = drainage_events(i,3); 
    % get image number of end of drainage, which is in column 3 of the
    % array 
    
    contem_lk = drainage_events(i,1); 
    % get lake number from column 1 of the array 
    
    contem_st_area = drainage_events(i,4); 
    % get SGL area at start of drainage 
    
    critical_refill = contem_st_area * (refill_amt_to_allow / 100)...
        * post_drain_cloudy_days; 
    % work out the critical refill threshold that is allowed within the
    % script - note that it takes into account the start area prior to
    % draiange, the number of post-drainage observations and the amount of
    % refilling that is allowed; the reason for multiplying the value by
    % the number of "post_drain_cloudy_days" is to allow it to refill to,
    % say, 18% on the first day, then remain at this value for the next 3
    % days (assuming that the post_drain_cloudy_days variable is set at
    % above, as otherwise this has no effect as it's just multiplied by 1
    % obviously), then the event would still be allowed. if you don't
    % include the "post_drain_cloudy_days" variable in calculation of the
    % "critical_refill" variable, then it would just sum all values in this
    % array and means that it would filter out events wrongly 
    
    subseq_day = contem_day + 1; 
    % get next day based upon the current day + 1 
    
    if subseq_day > image_number 
        drainage_events(i,11) = 5; % put a 5 here if so... 
    else
        
    post_drain_days = sampleSGLs(contem_lk,subseq_day:image_number); 
    % extract correct portion from this array based upon the current SGL
    % and the relevant number of images that remain afterwards 
    
    isnt_nan = sum(~isnan(post_drain_days)); 
    % find out how many subsequent cloud-free obs there are 
    drainage_events(i,11) = isnt_nan; % put this into column 11 
    
    post_drain_days(isnan(post_drain_days)) = []; % get rid of NaNs here
    % since we need cloud-free (or, here, no-data-free) observations 
    
    end
    
    if isnt_nan > (post_drain_cloudy_days - 1) % as long as the number of
        % cloud-free days is >number allowed minus 1, then: 
    
    post_drain_cloud_days = post_drain_days(1,1:post_drain_cloudy_days); 
    % extract the correct number from the array based upon the
    % post_drain_cloudy_days variable, to ensure we are not taking more
    % data than we actually want... 
    
    % this is a test based upon comparing two strings against each other to
    % see if they match; if they do, the value equals 1, so in other words
    % if we are not allowing refilling at all, then do something elseif do
    % something else 
    if strcmp(allow_refill,'N') == 1 % if refilling not allowed... 
    
        if sum(post_drain_cloud_days) == 0 
            % must all be 0s (no refilling) 
            
            drainage_events(i,12) = 1; % put a 1 in column 12 if 
            % condition met 
            
        else % if not met 
            
            drainage_events(i,12) = 0; % put a 0 in this column
            % (which I use to filter out the events later)  
            
        end
        
    elseif strcmp(allow_refill,'Y') == 1 % if refilling allowed...
        
        if sum(post_drain_cloud_days) <= critical_refill 
            % value must be <= the critical_refill value defined in
            % the script above 
            
            drainage_events(i,12) = 1; % put a 1 in this column if 
            % condition met  
            
        else % if condition not met 
            
            drainage_events(i,12) = 0; % put a 0 in column to filter later
            
        end % end of critical refill test just above 
    end % end of refilling allowed "if" test 
    end % end of "isnt_nan" test 
end % end of for-loop relating to number of events 

drainage_events_v1 = drainage_events; 
% make copy of this array prior to filtering 

drainage_events(drainage_events(:,12) == 0,:) = [];
% if the value in this column is 0, which indicates a "bad" event, as shown
% above, then we need to get rid of it, so make it a blank row by using []
% "bad" should probably "false positive" or "error of commission"

drainage_events(drainage_events(:,11) == 5,:) = [];

drainage_events_v2 = drainage_events; 
% make a new copy of the drainage_events variable prior to doing anything
% else with it (for checking purposes) 

%% Tidy up to remove duplicate identified events  

% Works by only taking the first row for a given SGL and not any of the
% other rows - acceptable 

unique_draining_SGLs = unique(drainage_events_v2(:,1)); 
% find out from the first column (i.e. SGL numbers) of the
% "drainage_events_v2" the number of unique entries - e.g. if column
% contains the following: [1 1 2 3 3 4 5 6 7 7], the output would be: 
% [1 2 3 4 5 6 7]. 

number_uniques = size(unique_draining_SGLs,1); 
% find out how many unique entries there are here and thus the number of
% unique rapidly draining SGLs 

drainage_events_v3 = zeros(1,12); 
% pre-allocating this variable to save time 

% do the following for the number of unique_draining_SGLs 
for i = 1:number_uniques
    
    current_unique = unique_draining_SGLs(i,1); 
    % get the number of the SGL 
     
    duplicate_events = (drainage_events_v2(:,1) == current_unique); 
    % make a new "logical" array where 1 means that the value matches the
    % "current_unique" in column 1 of the drainage_events_v2 array 
    
   	duplicate_drains = drainage_events_v2(duplicate_events(:,1) == 1,:); 
    % extract these duplicate drainage events based upon the new logical
    % array 
    
    drainage_events_v3(i,:) = duplicate_drains(1,:); 
    % put the details of this event into drainage events v3 array 
   
end 

%% Get only fast events (defined with threshold above) 

drainage_events_v4 = drainage_events_v3; 
% make a new copy of this array to avoid overrwriting 

drainage_events_v4(drainage_events_v4(:,8) > rapid_time_h,:) = []; 
% examine column 8, which contains the duration of the rapid drainage event
% in hours and compare this with the threshold allowed for rapid drainage -
% if the current value is more than this, then get rid of this by deleting
% the entire row 

%% Convert start & end days to Julian days 

number_events_new = size(drainage_events_v4,1); 
% find out the new length of the drainage_events_v4 array along the first
% dimension 

for i = 1:number_events_new % do the following for each event in the row
    
    curr_day_no_st = drainage_events_v4(i,2); 
    % current day of start is the second column and row is defined by the
    % loop index 
    
    curr_day_Julian_st = all_JD_nums(curr_day_no_st,1); 
    % this array called all_JD_nums will already be a part of the
    % workspace and will have the same order (although be in a column not
    % rows) because of the way the script works as the order of the images
    % in the big area_info array 
    % this means it's possible to just extract the relevant row from the
    % column and this represents the JD number of the start day of the
    % drainage event 
    
    drainage_events_v4(i,2) = curr_day_Julian_st; 
    % in this array, replace the start day with this JD of start as defined
    % above 
    
    % now do a nearly identical process as above but look at the third
    % column of the drainage_events_v4 array, which is the one containing
    % the end day number of the drainage event 
    curr_day_no_end = drainage_events_v4(i,3); 
    curr_day_Julian_end = all_JD_nums(curr_day_no_end,1); 
    drainage_events_v4(i,3) = curr_day_Julian_end;
    
end 

number_rapid_evts = size(drainage_events_v4,1); 

%% Make tidy array for rapidly draining SGLs & a column headings file 

rapid_drainage_events = drainage_events_v4(:,1:10); 
% extract every row but only extract rows 1 through to 10 because rows 11
% and 12 are redundant variables that were created above but which are not
% required for analysis 

% change directory and save output file 
cd(all_outputs); 
save('All_Rapid_Drainage_Events.mat','rapid_drainage_events'); 

% not writing to Excel file in this script version
% % find out if this array is empty or not
% emptiness = isempty(rapid_drainage_events); 
% 
% % if it's not empty, then write this to an Excel file 
% if emptiness == 0 
%     xlswrite('Rapid_Drainage_Events.xlsx',rapid_drainage_events);
% end

% make a cell containing the column headings (manually input) 
drainage_headings = {...                                            (Excel)
    'Lake number',...                                   % column 1    (A)
    'Start DOY',...                                     % column 2    (B)
    'End DOY',...                                       % column 3    (C)
    'Value at beginning (km2/m3)',...                   % column 4    (D)
    'Value at end (km2/m3)',...                         % column 5    (E)
    'Value change (km2/m3)',...                         % column 6    (F)
    'Value loss (%)',...                                % column 7    (G)
    'Event duration (h)',...                            % column 8    (H)
    'Mean DOY drainage date',...                        % column 9    (I)
    'Drainage date error (d)'};                         % column 10   (J)
 
% save these to a .mat file in all three of the new folders made 
save('Drainage_Headings.mat','drainage_headings'); 

cd(small_outputs); 
save('Drainage_Headings.mat','drainage_headings'); 

cd(big_outputs); 
save('Drainage_Headings.mat','drainage_headings'); 

% write out information to command line 
disp(['Details of rapidly draining lakes in "All_Rapid_Drainage_Events.mat"'...
    ' file and column headings in "Drainage_Headings.mat".']) 

%% Create arrays based upon SGL size (for S2-L8 vs. MODIS comparisons)

drainage_events_v5 = drainage_events_v4; 

if strcmp(drainage_type,'area') == 1 
    
% execute the following loop the number of times based upon the number of
% rapid events identified 
for i = 1:number_rapid_evts 
    
    % get the SGL number from the rapid_drainage_events array column 1 
    num_SGLake = rapid_drainage_events(i,1); 
    
    % use this number and go back to the original areas array and extract
    % all values from it 
    total_SGL_seasonal_areas = sampleSGLs(num_SGLake,:); 
    
    % find out the maximum area of this SGL over the course of the season -
    % to determine whether or not it fits into the category of SGLs that
    % would be detectable within MODIS imagery (i.e. at least 0.125 km2 in
    % size) 
    max_area = nanmax(total_SGL_seasonal_areas); 
    
    % need to get this maximum area back into pixels because I use the
    % variable of "tracking_area_big" in terms of pixels and not areas 
    max_pixels = max_area / (pix_size^2); 
   
    % if the maximum SGL area (in pixels) over the season is greater than
    % or equal to the "tracking_area_big" value (which is 139, since 139 *
    % 0.03 * 0.03 is just over 0.125 km2) 
    if max_pixels >= tracking_area_big 
        
        % then put a 1 into column 11 to indicate that it is a SGL with a
        % size that would be identifiable within MODIS imagery 
        drainage_events_v5(i,11) = 1; 
        
    else
        
        % if it not meeting this criterion for size, then the SGL must have
        % a maximum size that is smaller than that resolvable within MODIS
        % imagery, so put a 2 within this column, to be used for filtering
        % just below 
        drainage_events_v5(i,11) = 2; 
        
    end    
end 

% if the identification is based on volume thresholds, then do this... 
elseif strcmp(drainage_type,'volume') == 1 
       
% loading the lake areas for filtering drainage arrays based on size 
fname = [inputs 'lake_area_tracking.mat']; 
load(fname); 
lake_areas_filter = area_info_file; % specify variable accordingly 
final_row_num = size(lake_areas_filter,1); % identify final row 
lake_areas_filter(final_row_num,:) = []; % delete this row of the array

for i = 1:number_rapid_evts 
    
    num_SGLake = rapid_drainage_events(i,1); 
    all_season_areas = lake_areas_filter(num_SGLake,:); 
    max_area = nanmax(all_season_areas); 
    max_pixels = max_area / (pix_size^2); 
    if max_pixels >= tracking_area_big 
        drainage_events_v5(i,11) = 1; 
    else 
        drainage_events_v5(i,11) = 2; 
    end
end 
    
end % end of strcmp loop relating to type of lake drainage 

% make new arrays from the newly formed drainage_events_v5 array, which now
% contains information about maximum SGL size and its relation to the
% MODIS/Sentinel resolvable sizes 
small_rapid_drainage_events = drainage_events_v5; 
big_rapid_drainage_events = drainage_events_v5; 

% for the small array, we want to get rid of all examples which are equal
% to 1 since they are contained within the big events array 
small_rapid_drainage_events(small_rapid_drainage_events(:,11) == 1,:) = [];

% do the reverse for the big events array 
big_rapid_drainage_events(big_rapid_drainage_events(:,11) == 2,:) = [];

% get rid of this redundant column 11 for tidy array for analysis 
small_rapid_drainage_events = small_rapid_drainage_events(:,1:10); 
big_rapid_drainage_events = big_rapid_drainage_events(:,1:10); 

cd(small_outputs); 
save('Small_Rapid_Drainage_Events.mat','small_rapid_drainage_events');

% % find out if this array is empty or not
% emptiness2 = isempty(small_rapid_drainage_events); 
% 
% % if it's not empty, then write this to an Excel file 
% if emptiness2 == 0 
%     xlswrite('Small_Rapid_Drainage_Events.xlsx',...
%         small_rapid_drainage_events);
% end

cd(big_outputs); 
save('Big_Rapid_Drainage_Events.mat','big_rapid_drainage_events'); 

% % find out if this array is empty or not
% emptiness3 = isempty(big_rapid_drainage_events); 
% 
% % if it's not empty, then write this to an Excel file 
% if emptiness3 == 0 
%     xlswrite('Big_Rapid_Drainage_Events.xlsx',...
%         big_rapid_drainage_events);
% end

number_rpd_smalls = size(small_rapid_drainage_events,1); 
number_rpd_bigs = size(big_rapid_drainage_events,1); 

% display progress update 
disp('Acquired information about rapidly draining SGLs, so moving on to others.') 

secs = 1; pause(secs); % pause script for 1 seconds 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   End of rapidly draining SGLs section 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Remove rapidly draining SGLs from new area tracking array 

% find out all of the SGLs that drain rapidly by applying a unique function
% to the first row of this array 
rapid_SGLs = unique(rapid_drainage_events(:,1)); 

% find out the size of the number of rapid SGLs 
number_rapids = size(rapid_SGLs,1); 

% make a complete copy of the file that contains are information to allow
% the instances of rapidly draining SGLs to be removed from this array
% altogether 
copy_sampleSGLs = sampleSGLs; 

% for-loop to be repeated based upon the number of rapid SGLs 
for i = 1:number_rapids 
    
    % find out the number of the rapidly draining SGL from the first column
    curr_rapid_la = rapid_SGLs(i,1); 
    
    % remove this SGL from the new array that contains all other SGLs that
    % do not drain rapidly 
    copy_sampleSGLs(curr_rapid_la,:) = NaN;     
    
end 

%% Begin slow lake tracking process 

other_SGLs = copy_sampleSGLs; 

lakenum_other = size(other_SGLs,1); 

current_events_row = 0; % this needs setting at 0 outside of the for-loop
% as it increments by 1 within it once an event has been identified 

slw_drainage_events = zeros(1,10); 
% pre-allocate a base array of this size to speed up script 

for i = 1:lakenum_other 
    % do this for each of the lakes within the max extent 
    
    current_SGL_values = other_SGLs(i,:); 
    % get the current lake's values by extracting every column of current
    % row
    
    max_seasonal_SGL_area = nanmax(current_SGL_values); 
    % find out the maximum area of the SGL over the season (using the "max"
    % function may result in different output to "nanmax" here 
    
    critical_loss = max_seasonal_SGL_area * (slw_loss_reqd / 100); 
    % define the critical loss required to define a rapid SGL event based
    % upon the critical threshold inputted to the script and multiplying
    % the maximum seasonal area by this value 
    
    for p = 1:(image_number-1) % repeat this for-loop for every image but
        % don't do it for the last one as obviously there are no data
        % beyond it at that point to compare it with
        
        current_val = current_SGL_values(1,p); 
        % get the current SGL value by defining it using the "p" variable 
        
        next_values = current_SGL_values(1,(p + 1):end); 
        % extract everything from the current value + 1 to the end of the
        % array, meaning if the value is towards the start of season, there
        % are more data points against which this will be compared and so
        % the size of this array will vary over time as the for-loop
        % executes 
        
        number_nexts = size(next_values,2); 
        % find out how many "next_values" there are to compare the current
        % value against 
        
        for h = 1:number_nexts % do the following for each of the total
            % number of next values 
            
            next_val = next_values(1,h); 
            % get the current value being examined 
            
            if next_val <= (current_val - critical_loss) 
                % compare the current value against the next value and take
                % the difference between them; if the next_value is less
                % than the current_value minus the critical loss required
                % for a drainage event, then do the following 
                
                % write the following out to the command line - note that
                % it's only a preliminary event because these need tidying
                % up to avoid duplicates and so on (as I do below) 
%                  disp(['Found preliminary slow lake drainage event for ' ...
%                      'SGL number ' num2str(i) ...
%                      ' beginning on image number ' ...
%                      num2str(p) '.']); 
                
                % now put all of the relevant information relating to this
                % event into the drainage_events array: 
                current_events_row = current_events_row + 1; 
                % increment this variable by 1 to make sure current data
                % are not being overwritten 
                
               	slw_drainage_events(current_events_row,1) = i; 
                % put the SGL number into the first column 
                start_JDay = all_JD_nums(p,1); % get the JD number for 
                % this day (although this is not used here - the JD numbers
                % are added to this array later for various reasons and 
                % replace the current day numbers) 
                
                slw_drainage_events(current_events_row,2) = p;
                % put the current start day into column 2 
                
                % get the same information for the end day 
                end_day = h + p; 
                end_JDay = all_JD_nums(end_day,1); 
                slw_drainage_events(current_events_row,3) = end_day; 
                % put this information into row 3 
                
                slw_drainage_events(current_events_row,4) = current_val; 
                % find out the current value being examined on the start
                % day and put this into column 4 
                
                slw_drainage_events(current_events_row,5) = next_val; 
                % same for end area, put into column 5 
                
                change = current_val - next_val; 
                slw_drainage_events(current_events_row,6) = change; 
                % get the change in area and put this into column 6,
                % meaning that column 4 minus column 5 should equal column
                % 6 always
                
                percent_loss = 100 - (100 / current_val * next_val); 
                slw_drainage_events(current_events_row,7) = percent_loss; 
                % also get the percentage loss to ensure that events are
                % not being identified badly for whatever reason - put this
                % information into column 7 
                
                change_period = (end_JDay - start_JDay) * 24; 
                slw_drainage_events(current_events_row,8) = change_period;
                % work out period of time over which this event occurred by
                % taking end from start day and multiplying by 24 hours... 
                
                mean_drain_date = (start_JDay + end_JDay) / 2; 
                slw_drainage_events(current_events_row,9) = ...
                    mean_drain_date;
                % mean drainage date (I add JD number here at this stage as
                % this information is not used later) is calculated and put
                % into column 9 
                
                drain_date_error = (end_JDay - start_JDay) / 2; 
                slw_drainage_events(current_events_row,10) =...
                    drain_date_error;
                % error in drainage date is in column 10 
                
            end % end of if next_value < critical loss test          
        end % end of for each of the "next_values" 
    end % end of do it for each image in the dataset 
end % end of for each lake number in the dataset 
    
number_slw_events = size(slw_drainage_events,1);
% work out how many events there are in the array at the moment by taking
% the size along the first (i.e. rows) dimension 

slw_drainage_events_v0 = slw_drainage_events; 

%% Filter events where basin refills above critical threshold post-drainage

% repeat the following loop for each of the number_events identified in the
% array in total 
for i = 1:number_slw_events 
    
    contem_day = slw_drainage_events(i,3); 
    % get image number of end of drainage, which is in column 3 of the
    % array 
    
    contem_lk = slw_drainage_events(i,1); 
    % get lake number from column 1 of the array 
    
    contem_st_area = slw_drainage_events(i,4); 
    % get SGL area at start of drainage 
    
    rfill_slw_amt_allow = 100 - rfill_slw_amt_allow; 
    critical_refill = contem_st_area * (rfill_slw_amt_allow / 100)...
        * post_drain_cloudy_days; 
    % work out the critical refill threshold that is allowed within the
    % script - note that it takes into account the start area prior to
    % draiange, the number of post-drainage observations and the amount of
    % refilling that is allowed; the reason for multiplying the value by
    % the number of "post_drain_cloudy_days" is to allow it to refill to,
    % say, 18% on the first day, then remain at this value for the next 3
    % days (assuming that the post_drain_cloudy_days variable is set at
    % above, as otherwise this has no effect as it's just multiplied by 1
    % obviously), then the event would still be allowed. if you don't
    % include the "post_drain_cloudy_days" variable in calculation of the
    % "critical_refill" variable, then it would just sum all values in this
    % array and means that it would filter out events wrongly 
    
    subseq_day = contem_day + 1; 
    % get next day based upon the current day + 1 
    
    if subseq_day > image_number 
    
    slw_drainage_events(i,11) = 5;
    % put a 5 if this is the case as we need at least one cloud-free
    % observation of it, which we would not have in this case, so we can't
    % be sure of what's happened 
    
    else
    
    post_drain_days = other_SGLs(contem_lk,subseq_day:image_number); 
    % extract correct portion from this array based upon the current SGL
    % and the relevant number of images that remain afterwards 
    
    isnt_nan = sum(~isnan(post_drain_days)); 
    % find out how many subsequent cloud-free obs there are 
    slw_drainage_events(i,11) = isnt_nan; % put this into column 11 
    
    post_drain_days(isnan(post_drain_days)) = []; % get rid of NaNs here
    % since we need cloud-free (or, here, no-data-free) observations 
    
    end 
    
    if isnt_nan > (post_drain_cloudy_days - 1) % as long as the number of
        % cloud-free days is >number allowed minus 1, then: 
    
    post_drain_cloud_days = post_drain_days(1,1:post_drain_cloudy_days); 
    % extract the correct number from the array based upon the
    % post_drain_cloudy_days variable, to ensure we are not taking more
    % data than we actually want... 
    
    % this is a test based upon comparing two strings against each other to
    % see if they match; if they do, the value equals 1, so in other words
    % if we are not allowing refilling at all, then do something elseif do
    % something else 
    if strcmp(allow_slw_refill,'N') == 1 % if refilling not allowed... 
    
        if sum(post_drain_cloud_days) == 0 
            % must all be 0s (no refilling) 
            
            slw_drainage_events(i,12) = 1; % put a 1 in column 12 if 
            % condition met 
            
        else % if not met 
            
            slw_drainage_events(i,12) = 0; % put a 0 in this column
            % (which I use to filter out the events later)  
            
        end
        
    elseif strcmp(allow_slw_refill,'Y') == 1 % if refilling allowed...
        
        if sum(post_drain_cloud_days) < critical_refill 
            % value must b <= the critical_refill value defined in
            % the script above 
            
            slw_drainage_events(i,12) = 1; % put a 1 in this column if 
            % condition met  
            
        else % if condition not met 
            
            slw_drainage_events(i,12) = 0; % put a 0 in column to filter 
            % later
            
        end % end of critical refill test just above 
    end % end of refilling allowed "if" test 
    end % end of "isnt_nan" test 
end % end of for-loop relating to number of events 

slw_drainage_events_v1 = slw_drainage_events; 
% make copy of this array prior to filtering 

slw_drainage_events(slw_drainage_events(:,12) == 0,:) = [];
% if the value in this column is 0, which indicates a "bad" event, as shown
% above, then we need to get rid of it, so make it a blank row by using [] 

slw_drainage_events(slw_drainage_events(:,11) == 5,:) = [];
% also get rid of events where there is no post-drainage observation to
% confirm it 

slw_drainage_events_v2 = slw_drainage_events; 
% make a new copy fo the drainage_events variable prior to doing anything
% else with it (for checking purposes) 

%% Tidy up to remove duplicate identified events  

% Works by only taking the first row for a given SGL and not any of the
% other rows 

unique_draining_slw_SGLs = unique(slw_drainage_events_v2(:,1)); 
% find out from the first column (i.e. SGL numbers) of the
% "drainage_events_v2" the number of unique entries - e.g. if column
% contains the following: [1 1 2 3 3 4 5 6 7 7], the output would be: 
% [1 2 3 4 5 6 7]. 

number_uniques_slw = size(unique_draining_slw_SGLs,1); 
% find out how many unique entries there are here and thus the number of
% unique rapidly draining SGLs 

drainage_events_slw_v3 = zeros(1,12); 
% pre-allocating this variable to save time 

% do the following for the number of unique_draining_SGLs 
for i = 1:number_uniques_slw
    
    current_unique = unique_draining_slw_SGLs(i,1); 
    % get the number of the SGL 
     
    duplicate_events_slw = (slw_drainage_events_v2(:,1) == current_unique); 
    % make a new "logical" array where 1 means that the value matches the
    % "current_unique" in column 1 of the drainage_events_v2 array 
    
   	duplicate_drains_slw = ...
        slw_drainage_events_v2(duplicate_events_slw(:,1) == 1,:); 
    % extract these duplicate drainage events based upon the new logical
    % array 
    
    drainage_events_slw_v3(i,:) = duplicate_drains_slw(1,:); 
    % put the details of this event into drainage events v3 array 
   
end 

%% Get only slow events (defined with threshold above) 

drainage_events_slw_v4 = drainage_events_slw_v3; 
% make a new copy of this array to avoid overrwriting 

drainage_events_slw_v4(drainage_events_slw_v4(:,8) < slow_time_h,:) = []; 
% examine column 8, which contains the duration of the rapid drainage event
% in hours and compare this with the threshold allowed for rapid drainage -
% if the current value is more than this, then get rid of this by deleting
% the entire row 

%% Convert start & end days to Julian days 

number_events_slw_new = size(drainage_events_slw_v4,1); 
% find out the new length of the drainage_events_v4 array along the first
% dimension 

for i = 1:number_events_slw_new % do the following for each event in the row
    
    curr_day_no_st = drainage_events_slw_v4(i,2); 
    % current day of start is the second column and row is defined by the
    % loop index 
    
    curr_day_Julian_st = all_JD_nums(curr_day_no_st,1); 
    % this array called all_JD_nums will already be a part of the
    % workspace and will have the same order (although be in a column not
    % rows) because of the way the script works as the order of the images
    % in the big area_info array 
    % this means it's possible to just extract the relevant row from the
    % column and this represents the JD number of the start day of the
    % drainage event 
    
    drainage_events_slw_v4(i,2) = curr_day_Julian_st; 
    % in this array, replace the start day with this JD of start as defined
    % above 
    
    % now do an early identical process as above but look at the third
    % column of the drainage_events_v4 array, which is the one containing
    % the end day number of the drainage event 
    curr_day_no_end = drainage_events_slw_v4(i,3); 
    curr_day_Julian_end = all_JD_nums(curr_day_no_end,1); 
    drainage_events_slw_v4(i,3) = curr_day_Julian_end;
    
end 

number_slow_evts = size(drainage_events_slw_v4,1); 

%% Make tidy array for slowly draining SGLs & a column headings file 

slow_drainage_events = drainage_events_slw_v4(:,1:10); 
% extract every row but only extract rows 1 through to 10 because rows 11
% and 12 are redundant variables that were created above but which are not
% required for analysis 

% change directory and save output file 
cd(all_outputs); 
save('All_Slow_Drainage_Events.mat','slow_drainage_events'); 

% % find out if this array is empty or not
% emptiness4 = isempty(slow_drainage_events); 
% 
% % if it's not empty, then write this to an Excel file 
% if emptiness4 == 0 
%     xlswrite('Slow_Drainage_Events.xlsx',slow_drainage_events);
% end

% write out information to command line 
disp(['Details of slowly draining lakes in "Slow_Drainage_Events"'...
    ' file and column headings in "Drainage_Events_Headings".']) 

%% Create arrays based upon SGL size (for Sentinel vs. MODIS comparisons)

drainage_events_slw_v5 = drainage_events_slw_v4; 

if strcmp(drainage_type,'area') == 1 
    
% execute the following loop the number of times based upon the number of
% rapid events identified 
for i = 1:number_slow_evts 
    
    % get the SGL number from the rapid_drainage_events array column 1 
    num_SGLake = slow_drainage_events(i,1); 
    
    % use this number and go back to the original areas array and extract
    % all values from it 
    total_SGL_seasonal_areas = other_SGLs(num_SGLake,:); 
    
    % find out the maximum area of this SGL over the course of the season -
    % to determine whether or not it fits into the category of SGLs that
    % would be detectable within MODIS imagery (i.e. at least 0.125 km2 in
    % size) 
    max_area = nanmax(total_SGL_seasonal_areas); 
    
    % need to get this maximum area back into pixels because I use the
    % variable of "tracking_area_big" in terms of pixels and not areas 
    max_pixels = max_area / (pix_size^2); 
   
    % if the maximum SGL area (in pixels) over the season is greater than
    % or equal to the "tracking_area_big" value (which is 139, since 139 *
    % 0.03 * 0.03 is just over 0.125 km2) 
    if max_pixels >= tracking_area_big 
        
        % then put a 1 into column 11 to indicate that it is a SGL with a
        % size that would be identifiable within MODIS imagery 
        drainage_events_slw_v5(i,11) = 1; 
        
    else
        
        % if it not meeting this criterion for size, then the SGL must have
        % a maximum size that is smaller than that resolvable within MODIS
        % imagery, so put a 2 within this column, to be used for filtering
        % just below 
        drainage_events_slw_v5(i,11) = 2; 
        
    end    
end 

elseif strcmp(drainage_type,'volume') == 1 
       
% loading the lake areas for filtering drainage arrays based on size 
fname = [inputs 'lake_area_tracking.mat']; 
load(fname); 
lake_areas_filter = area_info_file; % specify variable accordingly 
final_row_num = size(lake_areas_filter,1); % identify final row 
lake_areas_filter(final_row_num,:) = []; % delete this row of the array

for i = 1:number_slow_evts
    
    num_SGLake = slow_drainage_events(i,1); 
    all_season_areas = lake_areas_filter(num_SGLake,:); 
    max_area = nanmax(all_season_areas); 
    max_pixels = max_area / (pix_size^2); 
    if max_pixels >= tracking_area_big 
        drainage_events_slw_v5(i,11) = 1; 
    else 
        drainage_events_slw_v5(i,11) = 2; 
    end
end 
    
end % end of strcmp loop relating to type of lake drainage 


% make new arrays from the newly formed drainage_events_v5 array, which now
% contains information about maximum SGL size and its relation to the
% MODIS/Sentinel resolvable sizes 
small_slw_drainage_events = drainage_events_slw_v5; 
big_slw_drainage_events = drainage_events_slw_v5; 

% for the small array, we want to get rid of all examples which are equal
% to 1 since they are contained within the big events array 
small_slw_drainage_events(small_slw_drainage_events(:,11) == 1,:) = [];

% do the reverse for the big events array 
big_slw_drainage_events(big_slw_drainage_events(:,11) == 2,:) = [];

% get rid of this redundant column 11 for tidy array for analysis 
small_slw_drainage_events = small_slw_drainage_events(:,1:10); 
big_slw_drainage_events = big_slw_drainage_events(:,1:10); 

cd(small_outputs); 
save('Small_Slow_Drainage_Events.mat','small_slw_drainage_events'); 

% find out if this array is empty or not
emptiness2 = isempty(small_slw_drainage_events); 
% 
% % if it's not empty, then write this to an Excel file 
% if emptiness2 == 0 
%     xlswrite('Small_Slow_Drainage_Events.xlsx',...
%         small_slw_drainage_events);
% end
cd(big_outputs); 
save('Big_Slow_Drainage_Events.mat','big_slw_drainage_events'); 

% find out if this array is empty or not
emptiness3 = isempty(big_slw_drainage_events); 

% % if it's not empty, then write this to an Excel file 
% if emptiness3 == 0 
%     xlswrite('Big_Slow_Drainage_Events.xlsx',...
%         big_slw_drainage_events);
% end

number_slow_smalls = size(small_slw_drainage_events,1); 
number_slow_bigs = size(big_slw_drainage_events,1); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   End of slowly draining SGLs section 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Make array to highlight potentially fast draining SGLs regarded as slow

% This is really probably not going to be used extensively as an output,
% but it creates another output based upon SGLs that drain "slowly" but
% which met some criteria for rapid drainage upon the day. Criteria that
% are met are "critical loss" threshold (of needing to lose more than 80%
% of the maximum seasonal area) and time period of change, but criterion
% which is not met is the next day not showing an increase in area of >=
% 20% since we then infer that the "middle day" of the three to be a
% misclassification and the SGL has not really drained rapidly, but it
% probably has drained slowly... 

% Not filtering out slowly drained SGLs based upon this but it might be
% % useful to have this variable stored for future 
% 
% % make a copy of this array 
% new_slows = slow_drainage_events; 
% final_possible = size(new_slows,1); % find out size along 1st dimension 
% 
% % repeat the number of loop so many times (probably until it's too late) 
% for i = 1:final_possible
%     
%     lost_area = new_slows(i,7); 
%     lost_time = new_slows(i,8); 
%     
%     if lost_area > rpd_loss_reqd && lost_time <= rapid_time_h 
%         
%         new_slows(i,11) = 1;
%         
%     else
%         
%         new_slows(i,11) = 0;
%     
%     end
% end 
% 
% slows_fasts_copy = new_slows; 
% slows_fasts_copy(slows_fasts_copy(:,11) == 0,:) = []; 
% 
% cd(all_outputs); 
% save('Potentially_Rapid_Lakes.mat','slows_fasts_copy'); 

disp(' ') 
disp('Finished getting information for slowly draining lakes.')

%% Final summary at command line 

disp(' ') 
disp('Final summary of lake drainage analysis:') 
disp(['Identified ' num2str(number_rapid_evts) ' rapid drainage events over all images.'])
disp(['Of this number, ' num2str(number_rpd_smalls) ' were for small SGLs (i.e. those unresolvable in MODIS).']) 
disp(['Of this number, ' num2str(number_rpd_bigs) ' were for big SGLs (i.e. those resolvable in MODIS).']) 
disp(' ') 
disp(['Identified ' num2str(number_slow_evts) ' slow drainage events or freeze-overs over all images.'])
disp(['Of this number, ' num2str(number_slow_smalls) ' were for small SGLs (i.e. those unresolvable in MODIS).']) 
disp(['Of this number, ' num2str(number_slow_bigs) ' were for big SGLs (i.e. those resolvable in MODIS).']) 
disp(' ') 
disp('Lake drainage analysis complete. Have fun!') 

cd(background_dir);

all_rapid_drainage_events = rapid_drainage_events; 
all_slow_drainage_events = slow_drainage_events; 
potentially_rapid_lakes = slows_fasts_copy; 

if strcmp(tidy_wkspace,'Y') == 1 
    
    clear area_info_file volume_info_file change change_period contem_day ...
        contem_lk contem_st_area copy_sampleSGLs critical_loss critical_refill ... 
        curr_date curr_day_Julian_end curr_day_Julian_st curr_day_no_end ... 
        curr_day_no_st curr_rapid_la current_events_row current_SGL_values ...
        current_unique current_val d dates day DOY drain_date_error ...
        drainage_events drainage_events_slw_v3 drainage_events_slw_v4 ...
        drainage_events_slw_v5 drainage_events_v1 drainage_events_v2 ...
        drainage_events_v3 drainage_events_v4 drainage_events_v5 drainage_type ...
        duplicate_drains duplicate_drains_slw duplicate_events duplicate_events_slw ...
        emptiness2 emptiness3 end_day end_JDay final_possible final_row_num fname ...
        h i image_number imagery inputs outputs isnt_nan lake_areas_filter lake_number ...
        lakenum_other lost_area lost_time max_area max_pixels max_seasonal_SGL_area ...
        mean_drain_date month new_slows next_val next_values num_SGLake number_events ...
        number_events_new number_events_slw_new number_nexts number_rapid_evts ...
        number_rapids number_rpd_bigs number_rpd_smalls number_slow_bigs ...
        number_slow_smalls number_slw_events number_uniques number_uniques_slw ...
        other_SGLs outputs p percent_loss pix_size post_drain_cloud_days ...
        post_drain_cloudy_days post_drain_days rapid_SGLs rapid_time_d ...
        rapid_time_h refill_amt_to_allow rfill_slw_amt_allow rpd_loss_reqd...
        sampleSGLs secs slow_time_d slow_time_h slw_drainage_events ...
        slw_drainage_events_v1 slw_draiange_events_v2 slw_loss_reqd small_outputs ...
        start_JDay subseq_day tracking_area_big background_dir ...
        unique_draining_SGLs unique_draining_slw_SGLs v v0 volume_info_file year ... 
        all_pitputs all_season_areas allow_refill allow_slw_refill big_outputs ...
        rapid_drainage_events slow_drainage_events slw_drainage_events_v0 ...
        slw_drainage_events_v2 all_outputs final_out_dir slows_fasts_copy ...
        number_slow_evts all_JD_nums 
    
end 
    
