Introduction
I was part of a group project that involved categorizing, grouping and labeling all the road intersections and interchanges of the county. I was responsible for developing the milestones of the project and the timeline for when each milestone needed to be started and completed. The first milestone of the project was to calculate the smallest angle formed by the roads of an intersection. Drawing on the ideas posted in this GIS Stack Exchange discussion, I created my own version of the workflow necessary to calculate the smallest angle of a road intersection.
The complexity of the overall project, and having several members on the group, made it beneficial to invest time developing a system of documentation to be used during the group project. These are some of the charts and graphics I created for our group.
GIS Work
The smallest angle for each intersection was calculated using the coordinates of the center of the intersection and the coordinates along each road, 30cm away from the center. A table containing all these coordinates was created from the outputs of many tools in ArcMap. These tool outputs were joined together, and the results saved as the Intersection_XY table.
The Intersection_XY table provides the data needed to find the angle between roads at an intersection. The process used involves first calculating the degree bearing that each road leaves the intersections. The difference between the degree bearing of any two roads is equal to angle in degrees between them.
The following Python code was used to calculate the bearing of each road as it leaves the center of the intersection. The math for the conversion of angle given by atan2 to bearing degree was derived by trial and error. The approach taken was to convert from degrees counter-clockwise (atan2’s output) to degrees clockwise (bearing) by adding or subtracting 180°. North starting at 0° required a rotation by 90°. The angle formed between any road of an intersection was then found by taking the difference between the bearings of those roads.
def get_bearing(xy1, xy2):
x1, y1 = xy1
x2, y2 = xy2
dx = x2 - x1
dy = y2 - y1
radian = math.atan2(dy, dx)
radian = radian - (math.pi / 2)
if (radian > 0):
degrees = 360 - (radian * 360) / (2 * math.pi)
else:
degrees = 360 - ((2 * math.pi + radian) * 360) / (2 * math.pi)
return degrees
bearing = get_bearing((center_x, center_y), (road_x, road_y))
Finally, I wrote a Python program that generates a graphical representation of any intersection using Matplotlib. It was written initially to make it easier to verify that the math used for calculating road bearing was accurate. The road intersection was recreated using a polar chart by graphing lines with the bearings that had been calculated. The polar chart was then visually compared with the same intersection viewed in ArcMap. This verification was done for a number of intersections until convinced that the process was working.
A text area in the lower left corner of the canvas allows a user to enter the intersection id in order to retrieve a different road intersection. Retrieving an intersection in this way was often easier than having to open ArcMap. The graphic could also be saved and used for a report.
