Solved! Go to Solution.
#assuming all features are single part!!! import arcpy lineFC = r"C:\temp\test.shp" dsc = arcpy.Describe(lineFC) shapeFieldName = dsc.ShapeFieldName featureNumber = 0 updateRows = arcpy.UpdateCursor(lineFC) for updateRow in updateRows: vertexList = [] shapeObj = updateRow.getValue(shapeFieldName) partObj = shapeObj.getPart(0) for pointObj in partObj: vertexList.append((pointObj.X, pointObj.Y)) vertexCount = len(vertexList) if vertexCount < 3: break #if only two verticies in the line, break the loop since there is no "angle" to calculate i = 0 # while i + 3 <= vertexCount: #that is until you run out of vertex triplicate pairs to count x1 = vertexList[0] y1 = vertexList[1] x2 = vertexList[i + 1][0] y2 = vertexList[i + 1][1] x3 = vertexList[i + 2][0] y3 = vertexList[i + 2][1] funkyStat = (x1 + y2) / 1.61803399 #blah blah whatever formula i = i + 1 updateRow.FUNKY_STAT = funkyStat updateRows.updateRow(updateRow) del updateRow, updateRows
#assuming all features are single part!!!
import arcpy
fcVertexList = []
lineFC = r"C:\temp\test.shp"
dsc = arcpy.Describe(lineFC)
shapeFieldName = dsc.ShapeFieldName
featureNumber = 0
searchRows = arcpy.SearchCursor(lineFC)
for searchRow in searchRows:
featureNumber = featureNumber + 1
fcVertexList.append([])
shapeObj = searchRow.getValue(shapeFieldName)
partObj = shapeObj.getPart(0)
for pointObj in partObj:
fcVertexList[featureNumber - 1].append([pointObj.X, pointObj.Y])
#assuming all features are single part!!! import arcpy lineFC = r"C:\temp\test.shp" dsc = arcpy.Describe(lineFC) shapeFieldName = dsc.ShapeFieldName featureNumber = 0 updateRows = arcpy.UpdateCursor(lineFC) for updateRow in updateRows: vertexList = [] shapeObj = updateRow.getValue(shapeFieldName) partObj = shapeObj.getPart(0) for pointObj in partObj: vertexList.append((pointObj.X, pointObj.Y)) vertexCount = len(vertexList) if vertexCount < 3: break #if only two verticies in the line, break the loop since there is no "angle" to calculate i = 0 # while i + 3 <= vertexCount: #that is until you run out of vertex triplicate pairs to count x1 = vertexList[0] y1 = vertexList[1] x2 = vertexList[i + 1][0] y2 = vertexList[i + 1][1] x3 = vertexList[i + 2][0] y3 = vertexList[i + 2][1] funkyStat = (x1 + y2) / 1.61803399 #blah blah whatever formula i = i + 1 updateRow.FUNKY_STAT = funkyStat updateRows.updateRow(updateRow) del updateRow, updateRows
Sure you can do it in the cursor... I guess an update cursor so you could then populate your angle statistic field value, right?
It'd look something like this then:
...
#
# 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
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!
field1 = "FIELD1"
field2 = "FIELD2"
updateRows = arcpy.UpdateCursor(fc)
for updateRow in updateRows:
field1Value = updateRow.getValue(field1)
updateRow.setValue(field2, field1Value / 2)
updateRows.updateRow(updateRow)
del updateRow, updateRows