This ntbk otlins a prcss that crtes mppngs betwn sgmnts f tw different GIS models depicting Colorado's roadway network. These models contain a variety of data, but none contain all the data needed to validate our statewide model. Given that all of these models represent the same network, mashing-up this data should be easy, right?

This exercise leverages ArcGIS's python SDK, and before analysis can begin various workspace paramters must be defined. The cell below define the working directories, the original shapefiles representing the networks and the spatial reference that will be used when outputting feature classes.


In [13]:
import arcpy as ap
import os


# create variables for directories that will be used throughout this notebook
local_directory = r'\\CDOTDBD01\Shared\Conflating_Networks'
directory_containing_original_shapefiles = local_directory + '\\network_shapefiles\Test'
local_geodatabase = local_directory + '\\working_files.gdb'


# define the workspace as the location of the original shapefiles and create a list of the shapefiles
ap.env.workspace = directory_containing_original_shapefiles
network_shapefiles = [x[:-4] for x in ap.ListFeatureClasses()]
ap.env.workspace = local_directory

# if the geodatabase doesn't exist, create it
if ap.Exists(local_geodatabase) == False:
    ap.CreateFileGDB_management(local_directory, 'working_files.gdb')
    
# define spatial reference as NAD_1983_2011_UTM_Zone_13N
spatial_reference = ap.SpatialReference(102382)
ap.env.outputCoordinateSystem = spatial_reference

ap.env.workspace = local_geodatabase

Mapping between links that start and/or stop in-between intersections is especially tricky. This arbitrary segmentation of networks can be observed below. This will be resolved by using the function defined below, which dissolves the entire network, then uses arcpy's feature to line tool to break lines wherever they intersect (i.e. intersections and overpasses). To run this function on our network shapefiles, uncomment the last two lines of the cell below.

Map 3. - Blue arrows represent the arbitrary segmentation of the TMC network, orange lines represent the DTD modeling network

In [2]:
def remove_arbitrary_segmentation_from(a_network_shapefile):

    # this function checks for a feature class, with the networks name appended with "_dissolved_to_line", if it exists the 
    # nothing happens, if it doesn't exist the original feature class is dissolved into one feature and then that feature is 
    # broken apart at each intersection of lines
    
    network_name = ap.Describe(directory_containing_original_shapefiles + '//' + a_network_shapefile + '.shp').baseName

    if ap.Exists('%s_dissolved_to_line_far' % network_name) == False:       
        ap.Dissolve_management(directory_containing_original_shapefiles + '//%s.shp' % network_name, 
                                   'in_memory\\%s_dissolved'% network_name,
                                   '',
                                   '',
                                   'MULTI_PART',
                                   'DISSOLVE_LINES')

        ap.FeatureToLine_management('in_memory\\%s_dissolved' % network_name,
                                            '%s_dissolved_to_line_far' % network_name)

        ap.Delete_management('in_memory\\%s_dissolved' % network_name)
Blw ar th to netwos ae bng conflated in this notebook. The Traffic Messaging Channel network (TMC) to the left, and CDOT's modeling network (DTD), to the right. Notice that only a fraction of the DTD modeling network is covered by the TMC network. The function in the next cell will export all links from each network that are within 20 meters of the other network. To run this function, uncmnt the lst three lns of the code cell blck blw.

Map 1. - The TMC network

Map 2. - DTD's modeling network


In [3]:
def create_feature_class_of_links_within_20m(network_1, network_2):
    
    # this function selects features from network_1 that are within 20 meters of network_2 and vice-versa
    
    for each_network in [network_1, network_2]:
        network_name = ap.Describe(each_network).baseName 
        if ap.Exists('%s_lyr' % network_name) == True:
            ap.Delete_management('%s_lyr' % ap.Describe(each_network).baseName)
        
        ap.MakeFeatureLayer_management(each_network, '%s_lyr' % network_name)
    
    network_1, network_2 = [ap.Describe(network_1).baseName,  network_name]
    
    ap.SelectLayerByLocation_management ('%s_lyr' % network_1, 'WITHIN_A_DISTANCE', '%s_lyr' % network_2, '20 Meters')
    ap.SelectLayerByLocation_management ('%s_lyr' % network_2, 'WITHIN_A_DISTANCE', '%s_lyr' % network_1, '20 Meters')
    
    for each_network in [network_1, network_2]:
        #if ap.Exists('%s' % each_network) == True:
            #ap.Delete_management(each_network)
        ap.CopyFeatures_management ('%s_lyr' % each_network, '%s' % each_network.replace('_far', ''))
        
        ap.Delete_management('%s_lyr' % each_network)
        ap.Delete_management(each_network)
Selecting features from a networks that are within 20 meters of the other network is an important step, but still includes lines that cross the other network. These lines are easily identified using arcpy's feature to vertices tool. The function used to create the black X identifying the dangling lines in Map 4 is created in the next cell.

Map 3. - DTD's modeling network lines that are within 20 meters of the TMC network

Map 4. - DTD's modeling network lines that are within 20 meters of the TMC network with dangles identified with Black X's

The functions in the next cell will use the feature class created in the cell above as an input and then create a feature class with dangling lines removed. The primary arcpy tool used in these functions is the feature vertices to points.

The first function below creates a feature class of dangling points from the feature class of links within 20 meters of the other network. Most of these 'new' dangling points should be removed, however some are not actually new dangles and represent the actual ends of the original network, in this case the dangle should not be removed. The second function, creates a feature class of the original dangling points. Finally, the third function compares the two to identify, then remove the new dangling lines.

To run this function on our network shapefiles, uncomment the last two lines of the cell below.

Next, intersection points will be created. Again, uncomment the last two lines to run this function.


In [4]:
def create_crossing_pts_from(a_network):
    network_name = ap.Describe(a_network).baseName
    
    for each_field in [a.name for a in ap.ListFields(a_network)]:
        try:
             ap.DeleteField_management(a_network,
                                      each_field)
        except:
            pass

    ap.FeatureVerticesToPoints_management ('%s' %  network_name, 
                                           '%s_pts' % network_name,
                                           'BOTH_ENDS')
    
    ap.DeleteIdentical_management('%s_pts' % network_name, 
                                  'Shape')

geomtric information from links will be created in the following cell, these attributes will later be transferred to intersections. The conflation will actually occur between intersections, rather than links.


In [5]:
import cmath, math

def the_coordinates_of(point_object):
    # function to return coordinates from arc point object
    x,y = round(point_object.X,5), round(point_object.Y,5)
    
    return x,y


def angle_between(point_a_coordinates, point_b_coordinates): 
    # function to calculate the angle between two sets of coordinatesdef angle_between(first_point_coordinates, last_point_coordinates):
    x1, y1 = float(point_a_coordinates[0]), float(point_a_coordinates[1])
    x2, y2 = float(point_b_coordinates[0]), float(point_b_coordinates[1])

    polarcoords = cmath.polar(complex(x2 - x1, y2 - y1))

    if math.degrees(polarcoords[1]) <= 90:
        degrees = 90 - math.degrees(polarcoords[1])
    else:
        degrees = 450 - math.degrees(polarcoords[1])

    return round(degrees, 1)



def the_coordinates_of_point_x_meters_from_the(start_or_end, of_a_line, int_of_x):
    # function returns coordinates of a point a number of meters from the start or end of the line
    # the first paramter must be a string of either 'start' or 'end'
    # second parameter is an arc line geometry 
    # and third is an int of the number of meters from either end of the line
    
    x_meters_as_percentage_of_length =  float(int_of_x)/float(of_a_line.length)
    
    if start_or_end.lower() == 'start':
        the_point = of_a_line.positionAlongLine(x_meters_as_percentage_of_length, True).firstPoint
    
    elif start_or_end.lower() == 'end':
        the_point = of_a_line.positionAlongLine(1 - x_meters_as_percentage_of_length, True).firstPoint
    
    else:
        print('''first parameter must be string of either "start" or "end"''')

    return the_coordinates_of(the_point)

            
              
def write_angles_from_start_of_each_line_of(a_line_feature_class):
    network_name = ap.Describe(a_line_feature_class).basename
    with ap.da.UpdateCursor('%s' % network_name,
            ['SHAPE@', 
            'first_10m_angle',
            'start_to_end_angle']) as a_cursor_for_line_angles:
          
        for each_line in a_cursor_for_line_angles:
            the_line_geometrys = of_a_line = each_line[0]

            first_point = the_coordinates_of(the_line_geometrys.firstPoint)
            last_point = the_coordinates_of(the_line_geometrys.lastPoint)
            point_10m_from_start = the_coordinates_of_point_x_meters_from_the('start', of_a_line, 10)
            
            angle_of_lines_first_10m = angle_between(first_point, point_10m_from_start)
            start_to_end_angle = angle_between(first_point,last_point)

            each_line[1] = str(angle_of_lines_first_10m)
            each_line[2] = str(start_to_end_angle)
              
            a_cursor_for_line_angles.updateRow(each_line)


def write_angles_from_end_of_each_line_of(a_line_feature_class):
    network_name = ap.Describe(a_line_feature_class).basename
    with ap.da.UpdateCursor('%s' % network_name,
            ['SHAPE@',  
            'last_10m_angle',
            'end_to_start_angle']) as a_cursor_for_line_angles:

          
        for each_line in a_cursor_for_line_angles:
            the_line_geometrys = of_a_line = each_line[0]

            first_point = the_coordinates_of(the_line_geometrys.firstPoint)
            last_point = the_coordinates_of(the_line_geometrys.lastPoint)
            point_10m_from_end = the_coordinates_of_point_x_meters_from_the('end', of_a_line, 10)
            
            angle_of_lines_last_10m =  angle_between(last_point, point_10m_from_end)
            end_to_start_angle = angle_between(last_point, first_point)

            each_line [1] = str(angle_of_lines_last_10m)
            each_line[2] = str(end_to_start_angle)
              
            a_cursor_for_line_angles.updateRow(each_line)

In [6]:
def create_start_points(shapefile):
    network_name = n = ap.Describe(shapefile).baseName
    ap.FeatureVerticesToPoints_management ('%s' % n,
                                           '%s_start_pts' % n,
                                           'START')
    
    
def add_start_angle_fields_to(line_fc):
    for each_field in ['first_10m_angle', 'start_to_end_angle']:
        
        if each_field not in [a.name for a in ap.ListFields(line_fc)]:
            ap.AddField_management(line_fc,
                  '%s' % each_field,
                  'STRING')


def create_start_points_with_angles(fc_of_lines):
    add_start_angle_fields_to('%s' % fc_of_lines)
    write_angles_from_start_of_each_line_of('%s' % fc_of_lines)
    create_start_points('%s' % fc_of_lines)

    

def join_start_points_to_all_points(start_pts, all_pts):
    ap.SpatialJoin_analysis(target_features= all_pts,
                                 join_features= start_pts, 
                                 out_feature_class="%s_to_start_pts" % all_pts, 
                                 join_operation="JOIN_ONE_TO_MANY", 
                                 join_type="KEEP_COMMON", 
                                 match_option="INTERSECT", 
                                 search_radius="1 Centimeters", 
                                 distance_field_name="")
    
    ap.Delete_management(start_pts)

In [7]:
def create_end_points(shapefile):
    network_name = n = ap.Describe(shapefile).baseName
    ap.FeatureVerticesToPoints_management ('%s' % n,
                                           '%s_end_pts' % n,
                                           'END')
    
    for fields in ['first_10m_angle', 'start_to_end_angle']:
        ap.DeleteField_management('%s_end_pts' % shapefile, 
                                  fields)
        
        
def add_end_angle_fields_to(line_fc):
    for each_field in [ 'last_10m_angle', 'end_to_start_angle']:
        
        if each_field not in [a.name for a in ap.ListFields(line_fc)]:
            ap.AddField_management(line_fc,
                  '%s' % each_field,
                  'STRING')

            
def create_end_points_with_angles(fc_of_lines):
    add_end_angle_fields_to(fc_of_lines)
    write_angles_from_end_of_each_line_of('%s' % fc_of_lines)
    create_end_points('%s' % fc_of_lines)        



def join_end_points_to_all_points(end_pts, all_pts):
    ap.SpatialJoin_analysis(target_features= all_pts,  
                            join_features=end_pts, 
                            out_feature_class="%s_to_end_pts" % all_pts, 
                            join_operation="JOIN_ONE_TO_MANY", 
                            join_type="KEEP_COMMON",
                            match_option="INTERSECT", 
                            search_radius="1 Centimeters", 
                            distance_field_name="")
    
    ap.Delete_management(end_pts)

In [8]:
from collections import defaultdict

def all_angles_to_dictionary_of_angles_with_line_id(all_pts):

    angle_dictionary = defaultdict(list)
    
    pts_joined_to_start = ('%s_to_start_pts' % all_pts, ["TARGET_FID", "first_10m_angle", "start_to_end_angle", "ORIG_FID_1"])
    pts_joined_to_end = ('%s_to_end_pts' % all_pts, ["TARGET_FID", "last_10m_angle", "end_to_start_angle", "ORIG_FID_1"])
      
    for each_join in [ pts_joined_to_start, pts_joined_to_end]:
        with ap.da.SearchCursor(each_join[0], each_join[1]) as cursor_to_find_angles:

            for row in cursor_to_find_angles:
                
                point_id = row[0]
                first_10m_angle = row[1]
                start_to_end_angle = row[2]
                line_id = row[3]
                
                if first_10m_angle != '<Null>':
                    try:
                        start_angle = float(first_10m_angle.strip())
                        line_angle = float(start_to_end_angle.strip())
                        angle_dictionary[row[0]].append((start_angle,line_angle, line_id))
                            
                    except:
                        continue
            
            ap.Delete_management(each_join[0])

    return angle_dictionary

def all_angles_to_dictionary_of_angles_without_line_id(all_pts):

    angle_dictionary = defaultdict(list)
    
    pts_joined_to_start = ('%s_to_start_pts' % all_pts, ["TARGET_FID", "first_10m_angle", "start_to_end_angle"])
    pts_joined_to_end = ('%s_to_end_pts' % all_pts, ["TARGET_FID", "last_10m_angle", "end_to_start_angle"])
      
    for each_join in [ pts_joined_to_start, pts_joined_to_end]:
        with ap.da.SearchCursor(each_join[0], each_join[1]) as cursor_to_find_angles:

            for row in cursor_to_find_angles:
                
                point_id = row[0]
                first_10m_angle = row[1]
                start_to_end_angle = row[2]
                
                if first_10m_angle != '<Null>':
                    try:
                        start_angle = float(first_10m_angle.strip())
                        line_angle = float(start_to_end_angle.strip())
                        angle_dictionary[row[0]].append((start_angle,line_angle))
                            
                    except:
                        continue
            
            ap.Delete_management(each_join[0])

    return angle_dictionary


def add_angles_from_dictionary_with_id_to_all_points(all_pts, angle_dictionary):
    for each_field in ["angle_to_10m_from_pt", "angle_to_end", "line_ids"]:
        ap.AddField_management(all_pts,'%s' % each_field,'STRING')
   

    with ap.da.UpdateCursor(all_pts, 
                            ["OBJECTID", "angle_to_10m_from_pt", "angle_to_end", "line_ids"]) as angle_cursor:
        
        for row in angle_cursor:
            point_id = row[0]
            
            angle_to_10m_from_pts = [x[0] for x in angle_dictionary[point_id]]
            angle_to_ends = [x[1] for x in angle_dictionary[point_id]]
            line_ids = [x[2] for x in angle_dictionary[point_id]]
            
            row[1] = str(angle_to_10m_from_pts).replace('[','').replace(']','')
            row[2] = str(angle_to_ends).replace('[','').replace(']','')
            row[3] = str(line_ids).replace('[','').replace(']','')
            
            angle_cursor.updateRow(row)
        
    del angle_dictionary, angle_cursor


def add_angles_from_dictionary_without_id_to_all_points(all_pts, angle_dictionary):
    for each_field in ["angle_to_10m_from_pt", "angle_to_end", "line_ids"]:
        ap.AddField_management(all_pts,'%s' % each_field,'STRING')
   

    with ap.da.UpdateCursor(all_pts, 
                            ["OBJECTID", "angle_to_10m_from_pt", "angle_to_end", "line_ids"]) as angle_cursor:
        
        for row in angle_cursor:
            point_id = row[0]
            
            angle_to_10m_from_pts = [x[0] for x in angle_dictionary[point_id]]
            angle_to_ends = [x[1] for x in angle_dictionary[point_id]]
           
            
            row[1] = str(angle_to_10m_from_pts).replace('[','').replace(']','')
            row[2] = str(angle_to_ends).replace('[','').replace(']','')

            
            angle_cursor.updateRow(row)
        
    del angle_dictionary, angle_cursor

In [9]:
import csv 

def locate_overpasses_from(original_network_fc):
    network_file = directory_containing_original_shapefiles + '//%s.shp' % original_network_fc
    network_name = ap.Describe(network_file).basename
    
    if ap.Exists('%s_to_line' % network_name) == False:
        ap.FeatureToLine_management(network_file,
                                    '%s_to_line' % network_name)
    
    for lyr in ['%s_to_line_lyr' % network_name ,'%s_dissolved_to_line_pts_lyr' % network_name]:
        if ap.Exists(lyr):
            ap.Delete_management(lyr)
        ap.MakeFeatureLayer_management(lyr[:-4], lyr)
    
    ap.SelectLayerByLocation_management ('%s_to_line_lyr' % network_name, 
                                         'WITHIN_A_DISTANCE', 
                                         '%s_dissolved_to_line_pts_lyr' % network_name,
                                         '1 Centimeters')
    
    list_of_overpasses = []
    with ap.da.SearchCursor('%s_to_line' % network_name, 
                           ['FID_%s' % network_name],
                            sql_clause = (None,'ORDER BY FID_%s' % network_name)) as overpass_cursor:

        last_feature_id = None
        for each_feature in overpass_cursor:
            feature_id = each_feature[0]

            if feature_id == last_feature_id:
                list_of_overpasses.append(last_feature_id)

            else:
                pass

            last_feature_id = feature_id
    
    
    ap.SelectLayerByAttribute_management('%s_to_line_lyr' % network_name, 
                                         'SUBSET_SELECTION',
                                         '"FID_%s" IN %s' % (network_name, tuple(list_of_overpasses)))
    
    
    ap.FeatureClassToFeatureClass_conversion ('%s_to_line_lyr' % network_name, 
                                              local_geodatabase, 
                                              '%s_overpass' % network_name)
    
    ap.Delete_management('%s_to_line_lyr' % network_name)
    ap.Delete_management('%s_dissolved_to_line_pts_lyr' % network_name)
    #ap.Delete_management('%s_to_line' % network_name)
    

def join_overpass_pts_to_dissolved_line_pts(dissolved_pts, overpass_pts):
    network_name = dissolved_pts.replace('_dissolved_to_line_pts','')
    ap.SpatialJoin_analysis(target_features= dissolved_pts,
                        join_features= overpass_pts,
                        out_feature_class="%s_pts_with_overpasses" % network_name, 
                        join_operation="JOIN_ONE_TO_MANY", 
                        join_type="KEEP_COMMON", 
                        match_option="INTERSECT", 
                        search_radius="1 Centimeters", 
                        distance_field_name="")
    
    pts_with_overpass = "%s_pts_with_overpasses" % each_network
    for each_field in  [a.name for a in ap.ListFields(pts_with_overpass)]:
       
        if each_field == 'angle_to_10m_from_pt_1':
            ap.AlterField_management (pts_with_overpass, 
                                      each_field, 
                                      "overpass_angle_to_10m", 
                                      "overpass_angle_to_10m_from_pt")

        if each_field =='angle_to_end_1':
            ap.AlterField_management (pts_with_overpass, 
                                      each_field,
                                      "overpass_angle_to_end", "overpass_angle_to_end")
            
        if each_field == 'angle_to_10m_from_pt':
            ap.AlterField_management (pts_with_overpass, 
                                      each_field, 
                                      "intersection_angle_to_10m", 
                                      "intersection_angle_to_10m_from_pt")

        if each_field =='angle_to_end':
            ap.AlterField_management (pts_with_overpass, 
                                      each_field,
                                      "intersection_angle_to_end", 
                                      "intersection_angle_to_end")
            
        if each_field =='line_ids':
            pass
        
        else:
            try:
                ap.DeleteField_management(pts_with_overpass, 
                                          each_field)
            except:
                pass
    
    #ap.Delete_management('%s_to_line' % network_name)
    
    
def update_intersections_from_overpasses(dissolved_pts_joined_to_overpass_pts):
    network_name = ap.Describe(dissolved_pts_joined_to_overpass_pts).baseName.replace('%s_pts_with_overpasses','')
    dict_for_new_points = {}
    
    with ap.da.SearchCursor(dissolved_pts_joined_to_overpass_pts, 
                            ['OBJECTID',
                            'intersection_angle_to_10m', 
                             'intersection_angle_to_end', 
                             'overpass_angle_to_10m',
                             'overpass_angle_to_end']) as intersection_cursor:
        
        for row in intersection_cursor:
            point_id = row[0]
            
            intersection_start_angles = row[1].split(',')
            intersection_line_angles = row[2].split(',')
            overpass_start_angles = row[3].split(',')
            overpass_line_angles = row[4].split(',')
            
            final_intersection_start_angles = [x for x in intersection_start_angles if x not in overpass_start_angles]
            final_intersection_line_angles = [x for x in intersection_line_angles if x not in overpass_line_angles]
            
            dict_for_new_points[point_id] = (final_intersection_start_angles, 
                                             final_intersection_line_angles,
                                             overpass_start_angles, 
                                             overpass_line_angles,)
            
    with open(r'\\CDOTDBD01\Shared\Conflating_Networks\%s_points_with_overpasses.csv' % network_name, 'wb') as csvfile:
        writer = csv.writer(csvfile, delimiter=',')
        writer.writerow(['point_id', 
                             'intersection_start_angles', 
                             'intersection_line_angles', 
                             'overpass_start_angles', 
                             'overpass_line_angles'])
            
        for point_id, angles in dict_for_new_points.items():
            i_start_angle, i_line_angle, o_start_angle, o_line_angle = angles
            
            writer.writerow([point_id, i_start_angle, i_line_angle, o_start_angle, o_line_angle])
    
    
    ap.MakeTableView_management(r'\\CDOTDBD01\Shared\Conflating_Networks\%s_points_with_overpasses.csv' % network_name, 'viewtemp')
    ap.TableToTable_conversion('viewtemp', local_geodatabase, '%s_angles_for_point' % network_name)
    ap.Delete_management('viewtemp')
    
    #ap.Delete_management('%s_overpass' % each_network)

In [10]:
def create_original_start_and_end_points_from(original_fc):
    name = ap.Describe(original_fc).baseName
    
    ap.MakeFeatureLayer_management(original_fc, '%s_lyr' % name)
    ap.MakeFeatureLayer_management('%s_dissolved_to_line_pts' % name, '%s_dissolved_to_line_pts_lyr' % name)
    
    ap.SelectLayerByLocation_management('%s_lyr' % name, 'WITHIN_A_DISTANCE', '%s_dissolved_to_line_pts_lyr' % name, '1 Centimeters')
    
    ap.FeatureToLine_management('%s_lyr' % name,
                                '%s_to_line' % name)
    
    
    create_start_points_with_angles('%s_to_line' % name)
    create_end_points_with_angles('%s_to_line' % name)
    
    for layers in ['%s_dissolved_to_line_pts_lyr', '%s_dissolved_to_line_pts_lyr']:
        ap.Delete_management(layers % name)

In [17]:
def overpasses_angles_to_dictionary(overpasses_pts):
    angles = defaultdict(list)
    
    for a in ['start', 'end']:
        if a == 'start':
            b = 'first'
            c = 'start_to_end_angle'
        else:
            b = 'last'
            c = 'end_to_start_angle'
        
        fields = ['TARGET_FID', '%s_10m_angle' % b, c]
        
        with ap.da.SearchCursor('%s_to_%s_pts' % (overpass_pts, a), 
                                fields) as cursor:
            for row in cursor:
                angles[row[0]].append((float(row[1]), float(row[2])))
                
                
    return angles

def final_intersections_from(overpass_angle_dictionary, to_dissolved_pts):
    overpass_pts = overpass_angle_dictionary.keys()
    dict_for_intersection_angles = {}
    fields = ['OBJECTID', 'angle_to_10m_from_pt', 'angle_to_end']
    
    with ap.da.SearchCursor(to_dissolved_pts, 
                                fields) as cursor:
            for row in cursor:
                if row[0] not in overpass_angle_dictionary.keys():
                    continue
                else:
                    overpass_start_angles = [x[0] for x in overpass_angle_dictionary[row[0]]]
                    overpass_line_angles =  [x[1] for x in overpass_angle_dictionary[row[0]]]
                    intersection_start_angles = row[1].split(',')
                    intersection_line_angles = row[2].split(',')
                    
                    intersection_start_angles_without_overpasses = [x for x in overpass_start_angles if x not in intersection_start_angles]
                    intersection_line_angles_without_overpasses = [x for x in overpass_line_angles if x not in intersection_line_angles]
                    dict_for_intersection_angles[row[0]] = (intersection_start_angles_without_overpasses, intersection_line_angles_without_overpasses)
    
    
    with open(r'\\CDOTDBD01\Shared\Conflating_Networks\%s_points_with_overpasses1.csv' % 'DTD_2010', 'wb') as csvfile:
        writer = csv.writer(csvfile, delimiter=',')
        writer.writerow(['point_id', 
                        'intersection_start_angles', 
                        'intersection_line_angles', 
                        'overpass_start_angles', 
                        'overpass_line_angles'])
        
        for point_id, int_angles in dict_for_intersection_angles.items():
            i_start_angle, i_line_angle = int_angles
            
            op_start_angles = str([x[0] for x in overpass_angle_dictionary[point_id]])
            op_line_angles = str([x[1] for x in overpass_angle_dictionary[point_id]])
            
            
            writer.writerow([point_id, i_start_angle, i_line_angle, op_start_angles, op_line_angles])

In [14]:
for each_network in network_shapefiles:
    remove_arbitrary_segmentation_from(each_network)

of_one = '%s_dissolved_to_line_far' % network_shapefiles[0]
another = '%s_dissolved_to_line_far' % network_shapefiles[1]

create_feature_class_of_links_within_20m(of_one, another)

In [15]:
for each_network in network_shapefiles:
   
    create_crossing_pts_from('%s_dissolved_to_line' % each_network)
    create_original_start_and_end_points_from(directory_containing_original_shapefiles + '\\' + each_network + '.shp')
   
    join_start_points_to_all_points('%s_to_line_start_pts' % each_network, '%s_dissolved_to_line_pts' % each_network)
    join_end_points_to_all_points('%s_to_line_end_pts' % each_network, '%s_dissolved_to_line_pts' % each_network)
    
    #angle_dictionary = all_angles_to_dictionary_of_angles_without_line_id('%s_dissolved_to_line_pts' % each_network)
    #add_angles_from_dictionary_without_id_to_all_points('%s_dissolved_to_line_pts' % each_network, angle_dictionary)
    
    #create_start_points_with_angles('%s_dissolved_to_line' % each_network)
    #create_end_points_with_angles('%s_dissolved_to_line' % each_network)

In [ ]:
for each_network in network_shapefiles:
    
    overpasses = '%s_overpass' % each_network
    overpass_pts = '%s_overpass_pts' % each_network
    overpass_start_pts = '%s_overpass_start_pts' % each_network
    overpass_end_pts = '%s_overpass_end_pts' % each_network
    
    locate_overpasses_from(each_network)
    create_crossing_pts_from(overpasses)
   
    create_start_points_with_angles(overpasses)
    join_start_points_to_all_points(overpass_start_pts, overpass_pts)
   
    create_end_points_with_angles(overpasses)
    join_end_points_to_all_points(overpass_end_pts, overpass_pts)
   
    #angle_dictionary = all_angles_to_dictionary_of_angles_without_line_id(overpass_pts)
    #add_angles_from_dictionary_without_id_to_all_points(overpass_pts, angle_dictionary)
    
    #join_overpass_pts_to_dissolved_line_pts('%s_dissolved_to_line_pts' % each_network, overpass_pts)
    
    #update_intersections_from_overpasses("%s_pts_with_overpasses" % each_network)

In [ ]:


In [16]:
for each_network in network_shapefiles:
    
    join_start_points_to_all_points('%s_dissolved_to_line_start_pts' % each_network, '%s_dissolved_to_line_pts' % each_network)
    join_end_points_to_all_points('%s_dissolved_to_line_end_pts' % each_network, '%s_dissolved_to_line_pts' % each_network)
    
    angle_dictionary = all_angles_to_dictionary_of_angles_with_line_id('%s_dissolved_to_line_pts' % each_network)
    add_angles_from_dictionary_with_id_to_all_points('%s_dissolved_to_line_pts' % each_network, angle_dictionary)

In [ ]:
a = [179.3,90.1,269.8,359.3]

b = [179.3,179.3,90.1,269.8, 369.3,359.3]

c =  [x for x in overpasses if x not in intersections]

print(c)

In [ ]:
a = define_intersections_from('DTD_2010_Links_with_directions_overpass_pt')

final_intersections_from(a, 'DTD_2010_Links_with_directions_dissolved_to_line_pts')

In [ ]:
for each_network in network_shapefiles:
    ap.CopyFeatures_management ('%s_pts_with_overpasses' % each_network, '%s_pts_with_overpasses_2' % each_network)
    with ap.da.UpdateCursor('%s_pts_with_overpasses_2' % each_network,
                           ['OBJECTID',
                             'Shape',
                             'intersection_angle_to_10m',
                             'intersection_angle_to_end',
                             'line_ids',
                             'overpass_angle_to_10m',
                             'overpass_angle_to_end']) as intersection_cursor:
        
        for row in intersection_cursor:

In [ ]:
[a.name for a in ap.ListFields('%s_pts_with_overpasses' % each_network)]
    
    
    #for each_fc in ['%s_overpass', '%s_overpass_pts', '%s_to_line']:
     #   ap.Delete_management(each_fc % each_network)

In [ ]:
for each_network in network_shapefiles:
    
    lines_from_pts = defaultdict(list)
    
    with ap.da.SearchCursor('%s_pts_with_overpasses' % each_network, 
                            ['OID@', 'Shape', 'line_ids', 'overpass_angle_to_10m']) as centerline_cursor:
        
        for row in centerline_cursor:
            point_id = row[0]
            point_xy = row[1]
            lines_from_pts[point_id].append(point_xy)
            
            
            list_of_line_ids = row[2].split(',')
            list_of_angles = [float(x.strip()) for x in row[3].split(',')]
            
            if len(list_of_angles) == len(list(set(list_of_angles))):
                for x in list_of_line_ids:
                
            else:

In [ ]:
if len(row[2]) == 1:
                continue
            else:
                point_id = row[0]
                angles_and_ids = [x.replace(' ','') for x in row[2:]]
                angles = angles_and_ids[1].split(',')

                all_angles_from_pt = [x.strip() for x in row[2].split(',')]
                angles_from_pt_without_duplicates = list(set(all_angles_from_pt))


                if len(angles) == len(list(set(angles))):
                    continue
                else:
                    one_way_angles = list(set([x for x in all_angles_from_pt if all_angles_from_pt.count(x) == 1]))
                    if len(one_way_angles) > 0:
                        for a in one_way_angles:
                            b = float(a)
                            lines_from_pts[point_id].append(b)
    
                    print(lines_from_pts[point_id])

In [ ]:
ap.SelectLayerByAttribute_management('%s_pts_with_overpasses' % each_network 
                                         'NEW_SELECTION',
                                         '"OID@" IN %s' % (tuple([x[0] for x in lines_from_pts])))

    with ap.da.SearchCursor('%s_pts_with_overpasses' % each_network, 
                            ['Shape', 'OID@']) as dual_carriageway_cursor:
        for row in dual_carriageway_cursor:
            point = row[0].firstPoint.X, row[0].firstPoint.Y
            for angles in lines_from_pts[row[0]]:

In [ ]:
for each_network in network_shapefiles:
    ap.MakeFeatureLayer_management('%s_to_line' % each_network, '%s_to_line_lyr' % each_network)
    ap.MakeFeatureLayer_management('%s_pts_with_overpasses' % each_network, '%s_dissolved_pts_lyr' % each_network)
    ap.SelectLayerByLocation_management ('%s_to_line_lyr' % each_network, 'WITHIN_A_DISTANCE', '%s_lyr' % each_network, '1 Centimeters')
    
    
    for a in ['%s_to_line_lyr','%s_dissolved_pts']:
        ap.Delete_management(a  % each_network)

In [ ]:


In [ ]:
for each_network in network_shapefiles:
    ids_of_centerlines = []
    with ap.da.SearchCursor('%s_overpass_pts' % each_network,
                            ['angle_to_10m_from_pt', 'line_ids']) as centerline_cursor:
        
        for row in centerline_cursor:
            list_of_angles = [x.strip() for x in row[0].split(',')]
            list_of_ids =  [x.strip() for x in row[1].split(',')]
            
            centerline = set()
            one_way = [x for x in list_of_angles if x not in centerline and not centerline.add(x)]
            
            for each_angle in list(centerline):
                centerline_id = list_of_ids[list_of_angles.index(each_angle)]
                ids_of_centerlines.append(centerline_id)
                
    
    #ap.SelectLayerByAttribute_management('%s_to_line_lyr' % network_name, 
     #                                    'NEW_SELECTION',
    #                                   '"FID_%s" IN %s' % (network_name, tuple(list_of_overpasses)))
    
    print((each_network, [int(x) for x in ids_of_centerlines]))

In [ ]:
def GetDifference_between(angle_1, angle_2):
    # Python modulus has same sign as divisor, which is positive here,
    # so no need to consider negative case
    r = (angle_2 - angle_1) % 360.0
    if r >= 180.0:
        r -= 360.0
    return (r)

len(list(set(!overpass_angle_to_end!.split(','))))