function [rs_cell,model_time] =...
    obs2rs_site_bounds(file_num,filepath,filename,site_index)
% Function will extract all information in regard to restoration site
% boundaries

for i = 1:file_num
    fnameloc(i) = {strcat(filepath, filename(i))};
end

polygon_27km = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'mass_Mo_chap3/abq_model_boundary_km0_27'],' ');

% Import observation file data
for j = 1:file_num
    if j == 1
        all_x_coords = (ncread(char(fnameloc{j}),'station_x_coordinate'))';
        all_y_coords = (ncread(char(fnameloc{j}),'station_y_coordinate'))';
        all_obs_ids_raw = (ncread(char(fnameloc{j}),'station_id'))';
        all_Vx = (ncread(char(fnameloc{j}),'x_velocity'))';
        all_Vy = (ncread(char(fnameloc{j}),'y_velocity'))';
        all_depth = (ncread(char(fnameloc{j}),'Waterdepth'))';
        all_wse = (ncread(char(fnameloc{j}),'waterlevel'))';
        time_stamp = ncread(char(fnameloc{j}),'time');
        max_time = max(time_stamp);
    elseif j > 1
        Vx_j = ncread(char(fnameloc{j}),'x_velocity');
        Vy_j = ncread(char(fnameloc{j}),'y_velocity');
        depth_j = ncread(char(fnameloc{j}),'Waterdepth');
        wse_j = ncread(char(fnameloc{j}),'waterlevel');
        time_stamp_j = ncread(char(fnameloc{j}),'time');
        
        all_Vx = [all_Vx; (Vx_j(:,time_stamp_j>max_time))'];
        all_Vy = [all_Vy; (Vy_j(:,time_stamp_j>max_time))'];
        all_depth = [all_depth; (depth_j(:,time_stamp_j>max_time))'];
        all_wse = [all_wse; (wse_j(:,time_stamp_j>max_time))'];
        time_stamp =...
            [time_stamp; time_stamp_j(time_stamp_j>max_time,:)];
        max_time = max(time_stamp);
        clear Vx_j Vy_j depth_j wse_j time_stamp_j
    end
end

%% Reduce all boundary points to those within 27 km polygon

in_polygon_yn = inpolygon(all_x_coords,all_y_coords,...
    polygon_27km(:,1),polygon_27km(:,2));

all_x_coords = all_x_coords(:,in_polygon_yn(:)==1);
all_y_coords = all_y_coords(:,in_polygon_yn(:)==1);
all_depth = all_depth(:,in_polygon_yn(:)==1);
all_obs_ids_raw = all_obs_ids_raw((in_polygon_yn(:)==1),:);
all_Vx = all_Vx(:,in_polygon_yn(:)==1);
all_Vy = all_Vy(:,in_polygon_yn(:)==1);
all_wse = all_wse(:,in_polygon_yn(:)==1);

%% Split all points into site boundaries

[~,num_sites] = size(site_index);
[~,obs_pts] = size(all_x_coords);
rs_cell = cell(num_sites,15);
rs_site_correction = [1 -1 -1 -1 1 1 -1 1 1 -1 1 -1 1 -1 -1 1 -1 1 1 1];

[time_steps,~] = size(time_stamp);

for s_num = 1:num_sites
    if s_num < 10;
        s_num_id = strcat('s',num2str(s_num),'_');
        rs_ids = cellstr(strcat(all_obs_ids_raw(:,1:3)));
        rs_yn = strncmp(rs_ids(:,1),s_num_id,3);
        rs_elems = sum(rs_yn(:));
        rs_i = 0;
        rs_elem_pt = zeros(rs_elems,2);
        rs_elem_Vx = zeros(rs_elems, time_steps);
        rs_elem_Vy = zeros(rs_elems, time_steps);
        rs_elem_depth = zeros(rs_elems, time_steps);
        
        for i = 1:obs_pts
            if rs_yn(i) == 1
                 rs_i = rs_i + 1;
                 rs_elem_pt(rs_i,1) = all_x_coords(1,i);
                 rs_elem_pt(rs_i,2) = all_y_coords(1,i);
                 rs_elem_Vx(rs_i,:) = all_Vx(:,i)';
                 rs_elem_Vy(rs_i,:) = all_Vy(:,i)';
                 rs_elem_depth(rs_i,:) = all_depth(:,i)';
            end
        end
    
    elseif s_num > 9
        s_num_id = strcat('s',num2str(s_num));
        rs_ids = cellstr(strcat(all_obs_ids_raw(:,1:3)));
        rs_yn = strncmp(rs_ids(:,1),s_num_id,3);
        rs_elems = sum(rs_yn(:));
        rs_i = 0;
        rs_elem_pt = zeros(rs_elems,2);
        rs_elem_Vx = zeros(rs_elems, time_steps);
        rs_elem_Vy = zeros(rs_elems, time_steps);
        rs_elem_depth = zeros(rs_elems, time_steps);
    
        for i = 1:obs_pts
            if rs_yn(i) == 1
                 rs_i = rs_i + 1;
                 rs_elem_pt(rs_i,1) = all_x_coords(1,i);
                 rs_elem_pt(rs_i,2) = all_y_coords(1,i);
                 rs_elem_Vx(rs_i,:) = all_Vx(:,i)';
                 rs_elem_Vy(rs_i,:) = all_Vy(:,i)';
                 rs_elem_depth(rs_i,:) = all_depth(:,i)';
            end
        end
    end
    
    rs_mid_point_xy = zeros(rs_elems-1,2);

    for p_i = 1:rs_elems-1
        rs_mid_point_xy(p_i,1) = (rs_elem_pt(p_i,1) +...
            rs_elem_pt(p_i+1,1)) / 2;
        rs_mid_point_xy(p_i,2) = (rs_elem_pt(p_i,2) +...
            rs_elem_pt(p_i+1,2)) / 2;
    end


    % Create matrices for boundary lengths based on element points

    % left_segs_pts
    % x_coord | y_coord | u | v | length | downstream_distance

    % Restoration site boundary calcs, boundary type = 1
    rs_normals = dflow_his_unit_vectors(1,rs_elem_pt);
    
    rs_normals(:,3) = rs_normals(:,3) * rs_site_correction(s_num);
    rs_normals(:,4) = rs_normals(:,4) * rs_site_correction(s_num);
    
    [rs_Q, rs_Mo_flux, rs_Mo_trans, rs_Mo_length, rs_area] =...
        dflow_mass_momentum_his_calc(1,rs_elem_Vx, rs_elem_Vy,...
        rs_elem_depth,time_steps,rs_normals);

    % Calculate distances upstream and downstream of each segment

    rs_distances = updownstream_distances(rs_normals,rs_elems);

    % Assign left and right information to cells
    % Cell 3 - element coordinates (NAD83 m)
    % Cell 4 - Vx (m/s)
    % Cell 5 - Vy (m/s)
    % Cell 6 - Depth (m)
    % Cell 7 - Discharge (m^3/s)
    % Cell 8 - Momentum flux per area (N/m^2)
    % Cell 9 - Momentum flux (N)
    % Cell 10 - Momentum length (m)
    % Cell 11 - Upstream and downstream distances (m)
    % Cell 12 - Cumulative downstream distance (m)
    % Cell 13 - Area of flux (m^2)
    
    rs_cell(s_num,1) = {s_num};
    rs_cell(s_num,2) = {rs_elems};
    rs_cell(s_num,3) = {rs_elem_pt};
    rs_cell(s_num,4) = {rs_elem_Vx};
    rs_cell(s_num,5) = {rs_elem_Vy};
    rs_cell(s_num,6) = {rs_elem_depth};
    rs_cell(s_num,7) = {rs_mid_point_xy};
    rs_cell(s_num,8) = {rs_Q};
    rs_cell(s_num,9) = {rs_Mo_flux};
    rs_cell(s_num,10) = {rs_Mo_trans};
    rs_cell(s_num,11) = {rs_Mo_length};
    rs_cell(s_num,12) = {rs_distances};
    rs_cell(s_num,13) = {rs_distances(:,1)};
    rs_cell(s_num,14) = {rs_area};
    rs_cell(s_num,15) = {rs_normals};

end
model_time = time_stamp;

end