% Extracts binary inundated information from last time step (typically from
% steady flow runs)

clear

%% Input files

folder_path = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012_histextent/run7023_f583cms/map_files/'];

polygon250m = load(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'attenuation_chap1/abq_valley_250m_polygons_0_25km']);

fnames = struct2cell(dir(strcat(folder_path,'*map.nc')))'; 
[fnum,~] = size(fnames); % number of files

for i = 1:fnum
    
    fname_i = strcat(folder_path,fnames{i,1});
    
    if i == 1
        id = double(ncread(fname_i,'FlowElemGlobalNr'));
        x = ncread(fname_i,'FlowElem_xcc');
        y = ncread(fname_i,'FlowElem_ycc');
        d = ncread(fname_i,'waterdepth');
        d_end = d(:,end);
        A = ncread(fname_i,'FlowElem_bac');
    else
        id = ([id; double(ncread(fname_i,'FlowElemGlobalNr'))]);
        x = ([x; ncread(fname_i,'FlowElem_xcc')]);
        y = ([y; ncread(fname_i,'FlowElem_ycc')]);
        d = ncread(fname_i,'waterdepth');
        d_end = ([d_end; d(:,end)]);
        A = ([A; ncread(fname_i,'FlowElem_bac')]);
    end
    
end

idxydA_raw = ([id x y d_end A]);
[~,unique_r,~] = unique(idxydA_raw(:,2:3),'rows','stable');
idxydA = idxydA_raw(unique_r,:);

inundated_yn = logical(idxydA(:,4)>0);

xs_poly_num = max(polygon250m(:,3));

area_sum = zeros(xs_poly_num,3);

for xs_num = 1:xs_poly_num
    
    in_poly_yn = inpolygon(idxydA(:,2),idxydA(:,3),...
        polygon250m(polygon250m(:,3)==xs_num,1),...
        polygon250m(polygon250m(:,3)==xs_num,2));
    
    area_sum(xs_num,1) = sum(idxydA(in_poly_yn==1,5));
    area_sum(xs_num,2) = sum(idxydA(in_poly_yn==1&inundated_yn==1,5));
    area_sum(xs_num,3) = sum(idxydA(in_poly_yn==1&inundated_yn==0,5));
    
end

