Sorry for late response - so much another work to do:( Ok, I understand, no problem with posting the code. Just to make it clear - it calculates "average curvature" of each line (angular change devided by length of line). I'm pretty sure there is a room for improvement in the script, I just did not have enough time to work on it...#
# This script calculates average "curvature" of roads
#
# It takes a line feature class (with pre-added fields "angleChng" and "AvCurv"!)
# then takes the fist line
# then takes first three verticies
# and calculates the angular change of line between those three verticies
# then it loops through all verticies of the line, calculates angular changes and cumulates them
# after going through whole line, it writes down the total angular change, calculates and writes down the average curvature
# and then this is repeated for all lines in the feature class
#
import arcpy, math, datetime, numpy
from arcpy import env
print "starting"
start = datetime.datetime.now() # for calculating time of process
lineFC = r"roads.shp" # input fc - with pre-added fields "angleChng" and "AvCurv"!
dsc = arcpy.Describe(lineFC)
shapeFieldName = dsc.ShapeFieldName
updateRows = arcpy.UpdateCursor(lineFC)
for updateRow in updateRows: # looping through lines
vertexList = []
shapeObj = updateRow.getValue(shapeFieldName)
partObj = shapeObj.getPart(0)
for pointObj in partObj:
vertexList.append((pointObj.X, pointObj.Y)) # getting coordinates of all line vertices
vertexCount = len(vertexList)
i = 0
angularChange=0
averageCurvature=0
while i + 3 <= vertexCount: # that is until you run out of vertex triplicate pairs to count
x1 = vertexList[0] # coordinates of three vertices
y1 = vertexList[1]
x2 = vertexList[i + 1][0]
y2 = vertexList[i + 1][1]
x3 = vertexList[i + 2][0]
y3 = vertexList[i + 2][1]
u1=x1-x2 # vectors (BA and BC !)
u2=y1-y2
v1=x3-x2
v2=y3-y2
cit=(u1*v1)+(u2*v2) # calculating the angle of vectors
jmen1=math.sqrt((u1*u1)+(u2*u2))
jmen2=math.sqrt((v1*v1)+(v2*v2))
jmen=jmen1*jmen2
cosfi=cit/jmen
fiRad=numpy.arccos(cosfi) # angle of vectors [rad]
fiGrad=fiRad*200/math.pi
fiGrad=200-fiGrad # angle change of one segment of line [grad]
angularChange=angularChange+fiGrad # total angle change of line
i = i + 1
averageCurvature=angularChange/shapeObj.length*1000 # average curvature of line = angle change in one kilometer
updateRow.angleChng=angularChange # adding angle change of line
updateRow.AvCurv=averageCurvature # adding average curvature of line
updateRows.updateRow(updateRow)
print "Done in ",datetime.datetime.now() - start, " seconds"
del updateRows, updateRow, lineFC
Btw.: Maybe a newbie question, but i didn't find any solution yet - Is it possible to add field into fc (with fieldName as a variable) and the to use something like updateRow.fieldName? (if using posted script, it is necessary to have "pre-added" fields - it would be much easier to add them in the beginning of the script...) Thanks!
