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