# Finding angle between vertices to estimate curvature in the same feature class using python script?

Background info

What I'm attempting to create is an estimation of curvature, by finding the angle between vertices. The vertices layer was created from using the "Split line at vertices tool" on a roads layer. From there I used calculate geometry to get xy data for the line.

The problem

Create a function which uses a line's start and endpoint xy data to determine the angle between it's neighboring line segment(s). What I'm struggling with is determining what how the lines connect. The angle can be determined using math.atan2(), but I have no idea how find the connecting segments.

x1- x coordinate for the start of the line segment

y1- y coordinate for the start of the line segment

x2 -x coordinate for the end of the line segment

y2- y coordinate for the end of the line segment

Basic Pseudo code

``````function angle(x1,y1,x2,y2)

Given x1,y1 in column[observed line segment]

Find x2,y2 in column[neighboring line segment] == x1, y1 in column[observed line segment]

If lines segments have matching coordinates find angle between lines.

return angle
``````
• can you add some more information about what tools you're using? Are you working with ArcMap or QGIS? feature class? shapefile? lines or polygons?
– crld
Jun 17, 2016 at 17:30
• If you have a Standard or Advanced ArcMap license, you could try the "Update COGO Attributes. This will provide the direction of the line. I think you need to Add COGO attributes first. Jun 17, 2016 at 17:30
• Just work with line vertices Jun 17, 2016 at 19:08

I don't think you need to do your initial line splitting. If you iterate through the points in the `SHAPE@JSON` attribute of the feature class, you can measure the sum of angles for each line and divide by the number of vertices.

``````import re
import math
import arcpy
from math import acos
from numpy.linalg import norm
import numpy
import json

def findangle(p1,p2,p3):

"""This function returns the angle p1p2p3 given tuple inputs of the coordinates."""

#The vectors from p2->p1 and p2->p3
v1 = [p1-p2, p1-p2]
v2 = [p3-p2, p3-p2]

#The sum of the dotproducts of each vector.
dotproduct = numpy.sum(p*q for p,q in zip(v1,v2))

#The multiplication of the magnitudes of each vector.
magnitude = float(norm(v1) * norm(v2))

#Rounding here to account for number of decimal places in rads.
if dotproduct != 0 and magnitude != 0.0:
try: x = dotproduct/magnitude
except: print dotproduct, magnitude

if x > 1.0 or x <-1.0:
x = round(x,1)
if x > 1.0 or x < -1.0:
x = float(int(x))
return math.degrees(math.acos(round((x), 4)))
else: return 0
``````

Below is an example of how to measure the angles of vertices in a polygon

``````    with arcpy.da.SearchCursor(polygons, polygonFields) as cursor:
for row in cursor:
count +=1
if row is not None:
angle_list = []
for ringIndex, ring in enumerate(ringlist):
For measuring the angles of a line, the difference would be the `['rings']` value. You should print the `SHAPE@JSON`string for a line and look at the key that has all the vertices in it. Also different: polygon 'rings' are essentially lines where the first and last points are the same. In the logic above, if the current point is the first point of the line, it measures the angle from the next point, current point, and the second to last point. For your purposes, you might want to change that to just use the last point, or skip it if that information isn't useful.