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

Question

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 

Answer 1

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[0]-p2[0], p1[1]-p2[1]]
    v2 = [p3[0]-p2[0], p3[1]-p2[1]]

    #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[1] is not None:
                ringlist = json.loads(row[1].JSON)["rings"]
                angle_list = []
                for ringIndex, ring in enumerate(ringlist):
                    vertexList = json.loads(row[1].JSON)["rings"][ringIndex]
                    for index, vertex in enumerate(vertexList):
                        [x,y] = vertex
                        ptID = index
                        ringID = ringIndex
                        if ptID == 0: angle = findangle(vertexList[len(vertexList)-2], vertexList[0], vertexList[1])
                        elif ptID == len(vertexList)-1: angle = findangle(vertexList[index-1], vertexList[index], vertexList[1])
                        else: angle = findangle(vertexList[index-1], vertexList[index], vertexList[index+1])
                        angle_list.append(angle)
                 angle_sum = sum(angle_list)
    del cursor

For measuring the angles of a line, the difference would be the ['rings'] value. You should print the SHAPE@JSONstring 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.

arcpy.Geometry (http://desktop.arcgis.com/en/arcmap/10.3/analyze/arcpy-classes/geometry.htm)