% Extracts information from DFlow to points associated with a line boundary
% digitized in GIS to represent a channel-floodplain boundary

clear

%% ***INPUT*** Load polygon and polyline vertices 

disp('Collecting polygon and polyline information...')

% Input polyline file name and location
polyline_floc = '/Users/Colin/Documents/MATLAB/DFlow/xs_pt_extraction/';
polyline_fname = 'xs_250m_split_vertices';
polyline_xyij = importdata(strcat(polyline_floc, polyline_fname),' ');

% Input polygon file name and location
exist_polygon = 1; % Does the polygon exist? 1 - yes; 0 - no
polygon_floc = '/Users/Colin/Documents/MATLAB/DFlow/xs_pt_extraction/';
polygon_fname = 'xs_250m_buffer_polygons';
polygon_fnameloc = strcat(polygon_floc, polygon_fname);

%% ***INPUT*** Load DFlow data using partmap2variable function (raw data) 
% Load D-Flow information file either already turned to text file or raw
% D-Flow output if desired

disp('Collecting DFlow data from partitions...')

% Input file location, partition-0000 file name, and # partitions
part_num = 32;
dflow_floc = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3034_s180pback/'];
dflow_beg = 'abq10a12bfs180pback_';
dflow_end = '_20160101_000000_map.nc';

all_elem_data = partmap2elements(part_num,dflow_floc,dflow_beg,dflow_end);
    
polygon_data = importdata(polygon_fnameloc,' ');
x_polygon = polygon_data(:,1);
y_polygon = polygon_data(:,2);

xs_num = max(polygon_data(:,3));

for xs_i = 1:xs_num

%% Use a polygon to cut down on number of data points
% (polygon nodes, element ID, element x & y)

    if exist_polygon == 1
        disp('Checking which elements are inside of polygon...')

        % Extract proper elements for each xs polygon
        
        x_polygon_i = x_polygon(polygon_data(:,3)==xs_i,:);
        y_polygon_i = y_polygon(polygon_data(:,3)==xs_i,:);

        elems_in_poly_check = inpolygon(...
            all_elem_data(:,3),all_elem_data(:,4),...
            x_polygon_i,y_polygon_i);
        elems_in_poly = all_elem_data(elems_in_poly_check(:,1)==1,:);

    elseif exist_polygon == 0
        elems_in_poly = all_elem_data;
    end

    %% Find bounding nodes of elements within polygon
    disp('Retrieving bounding nodes...')

    [element_nodes_id, element_nodes_x, element_nodes_y] = ...
        find_bounding_nodes(dflow_floc, dflow_beg, dflow_end,...
        elems_in_poly);
    
    for j = 1:3

    %% Using a boundary polyline determine if element is intersected
    % (polyline nodes, element ID, bounding node x & y)

        polyline_xy = polyline_xyij(polyline_xyij(:,3)==xs_i & ...
            polyline_xyij(:,4)==j,1:2);
    
        disp(['Determining which elements the boundary polyline'...
            'intersects...'])
        [line_yes_no] = line_elem_intersect(element_nodes_x,...
            element_nodes_y,polyline_xy);

        elems_bndry = elems_in_poly(line_yes_no(:,1)==1,:);
        elem_bndry_nodes_x = element_nodes_x(line_yes_no(:,1)==1,:);
        elem_bndry_nodes_y = element_nodes_y(line_yes_no(:,1)==1,:);

        %% Find northern-most element of intersected elements
        % (element ID, element x & y)
        disp('Finding most western element...')

        elem_x = elems_bndry(:,3);
        elem_y = elems_bndry(:,4);

        [min_elem_x,min_elem_index] = min(elem_x);

        %% Order elements from upstream to downstream 
        % Link elements using shared bounding nodes - connected elements
        % will share 2 nodes (element ID, element x & y, bounding node x &
        % y)

        disp('Ordering elements right to left...')

        [elems_ordered, nodes_x_order, nodes_y_order,elems_remain] =...
            xs_element_order(min_elem_index, elems_bndry,...
            elem_bndry_nodes_x, elem_bndry_nodes_y);

        disp('Ordered.')
        
        if xs_i == 1 && j == 1
            
            xs_elements(:,1:5) = elems_ordered;
            xs_elements(:,6) = xs_i*10 + j;
            
        else
            xs_elem_hold(:,1:5) = elems_ordered;
            xs_elem_hold(:,6) = xs_i*10 + j;
            xs_elements = [xs_elements; xs_elem_hold];
        end
        
        clear xs_elem_hold
        
    end
    
    fprintf('%% complete: %d\r',xs_i/xs_num*100)

end