This notebook outlines a process that creates mappings between segments of two 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.
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)
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)
In [ ]:
## Shelley's scratch
# Import arcpy module
import arcpy
# Local variables:
NAV_TEST2 = "NAV_TEST2"
DTD_TEST2_FeatureVerticesToP1 = "DTD_TEST2_FeatureVerticesToP1"
DTD_TEST2_FeatureVerticesToP1__2_ = DTD_TEST2_FeatureVerticesToP1
NAV_TEST2__2_ = "NAV_TEST2"
NAV_TEST2_FeatureVerticesToP2 = "C:\\Users\\broadways\\Documents\\ArcGIS\\Default.gdb\\NAV_TEST2_FeatureVerticesToP2"
NAV_TEST2_FeatureVerticesToP1 = "NAV_TEST2_FeatureVerticesToP1"
NAV_TEST2_FeatureVerticesToP1__2_ = NAV_TEST2_FeatureVerticesToP1
DTD_TEST2 = "DTD_TEST2"
DTD_TEST2_FeatureVerticesToP = "C:\\Users\\broadways\\Documents\\ArcGIS\\Default.gdb\\DTD_TEST2_FeatureVerticesToP"
# Process: Feature Vertices To Points (2)
arcpy.FeatureVerticesToPoints_management(NAV_TEST2__2_, NAV_TEST2_FeatureVerticesToP2, "MID")
# Process: Select Layer By Location
arcpy.SelectLayerByLocation_management(DTD_TEST2_FeatureVerticesToP1, "WITHIN_A_DISTANCE", NAV_TEST2_FeatureVerticesToP2, "20 Meters", "NEW_SELECTION", "NOT_INVERT")
# Process: Feature Vertices To Points
arcpy.FeatureVerticesToPoints_management(DTD_TEST2, DTD_TEST2_FeatureVerticesToP, "MID")
# Process: Select Layer By Location (2)
arcpy.SelectLayerByLocation_management(NAV_TEST2_FeatureVerticesToP1, "WITHIN_A_DISTANCE", DTD_TEST2_FeatureVerticesToP, "20 Meters", "NEW_SELECTION", "NOT_INVERT")
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(','))))