[GIS] Calculating angle of polyline at midpoint using Python/ArcPy

arcgis-desktoparcpymathematicspython

Trying to calculate the angle of a line at its mid point using python; this will be used to display flow direction on a water network.

I've developed this script, which works in terms of generating the correct output, however, the angle returned is incorrect. I am currenlty getting some negative values, where I would like to have outputs of 0 to 360, in either an Arithmetic or Geographic format.

#OperationAlligator.py
import arcpy, math, datetime, numpy
from arcpy import env

print ("starting")
start = datetime.datetime.now() # for calculating time of process

#setting the test sources
flowArrow = r"J:\PYTHON\Flow_Direction.gdb\test" # pre-added fields "angle"
pipes = r"J:\PYTHON\Flow_Direction.gdb\Pipes"

#UpdateCursor
rows = arcpy.UpdateCursor(pipes)
updateRows_FlowArrow = arcpy.UpdateCursor(flowArrow)


#housekeeping
arcpy.DeleteFeatures_management(flowArrow)
arcpy.env.overwriteOutput = True


#Creating the physical xy locations in memory
shapeName = arcpy.Describe(pipes).shapeFieldName
for row in rows:
    feat = row.getValue(shapeName)
    x = feat.positionAlongLine(0.51,True).firstPoint.X - feat.positionAlongLine(0.49,True).firstPoint.X
    y = feat.positionAlongLine(0.51,True).firstPoint.Y - feat.positionAlongLine(0.49,True).firstPoint.Y
    radian = math.atan2(y,x)

    degrees = radian * 180 / math.pi
    print degrees
    row.angle = degrees
    rows.updateRow(row)

#Generating the mid point flow arrow
    out_cursor = arcpy.InsertCursor(flowArrow)
    midpoint = feat.positionAlongLine(0.50,True).firstPoint
    outRow = out_cursor.newRow()
    outRow.Shape = midpoint
    outRow.angle = row.angle
    out_cursor.insertRow(outRow)

del rows, row, updateRows_FlowArrow, outRow, out_cursor
print "Done in ",datetime.datetime.now() - start, " seconds"

visual display of the midpoint angle calculation problem

Best Answer

Try replacing

radian = math.atan2((x1 - y1),(x2 - y2))

with

radian = math.atan2((y2 - y1),(x2 - x1))

You should also avoid using the ends of your line segments but rather use the points located just before and after your midpoint

radian = math.atan2((feat.positionAlongLine(0.51,True).firstPoint.Y - feat.positionAlongLine(0.49,True).firstPoint.Y),(feat.positionAlongLine(0.51,True).firstPoint.X - feat.positionAlongLine(0.49,True).firstPoint.X))

Related Solutions

[GIS] Creating point at midpoint of polyline with ArcPy

this should do it in arcgis...

import arcpy, math, datetime, numpy

print ("starting")
start = datetime.datetime.now() # for calculating time of process

#setting the containers
midpoint = #put your file to be populated in here, make sure it already exists
polyline = #put your polyline file in here

#housekeeping
arcpy.DeleteFeatures_management(midpoint)
arcpy.env.overwriteOutput = True

#generating the mid point
with arcpy.da.SearchCursor(polyline, "SHAPE@") as in_cursor, \
arcpy.da.InsertCursor(midpoint, "SHAPE@") as out_cursor:
for row in in_cursor:
    midpoint = row[0].positionAlongLine(0.50,True).firstPoint
    out_cursor.insertRow([midpoint])

#tidy up
del rows, row, updateRows_midpoint, outRow, out_cursor, midpoint, polyline
print "Done in ",datetime.datetime.now() - start, " seconds"

[GIS] Calculating angle between polyline segments and assign this value to point

Here's a little script I put together to read a line and write points at vertices:

import os, sys, arcpy, math
# read http://help.arcgis.com/En/Arcgisdesktop/10.0/Help/index.html#/Adding_a_script_tool/00150000001r000000/
InLines = sys.argv[1] # input:Feature Class
OutFile = sys.argv[2] # output:Shapefile

# using OS functions:
# assume OutFile is "C:\\this\\path\\file.shp"
# then os.path.dirname(OutFile) = C:\this\path
# and os.path.basename(OutFile) = file.shp
# http://help.arcgis.com/En/Arcgisdesktop/10.0/Help/index.html#//00170000002p000000
arcpy.CreateFeatureclass_management(os.path.dirname(OutFile),os.path.basename(OutFile),"POINT")
arcpy.AddField_management(OutFile,"ANGLE","DOUBLE") # Add the angle field or there's no point doing this..

# you might want to do something like this:
# SR = arcpy.SpatialReference(28356) # use EPSG/SRID code here to create a spatial reference
# arcpy.DefineProjection_management(OutFile,SR), but you can do that later..

# create an insert cursor to write the output and a search cursor to read the input
with arcpy.da.InsertCursor(OutFile,["SHAPE@XY","ANGLE"]) as ICur:
    with arcpy.da.SearchCursor(InLines,"SHAPE@") as SCur:
        for ThisFeature in SCur:
            # loop for every feature..
            # read the geometries like http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//002z0000001t000000
            Geom = ThisFeature[0] # should be a polyline but polygons will work too
            Parts = range(Geom.partCount) #http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//000v000000n2000000
            # polygons and polygons can be made of many parts
            # so it's important to read each part (it's usually only 1 part)
            for prt in Parts:
                ThisPart = Geom.getPart(prt) # ThisPart is a Geometry object
                PointCount = ThisPart.count  # http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Array/000v0000005r000000/
                qPointsRange = range(1,PointCount -1) # range from 2nd to 2nd last

                for point in qPointsRange:
                    ThisPoint = ThisPart.getObject(point) # ThisPoint is a Point object http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Point/000v000000mv000000/
                    X = ThisPoint.X # this vertex X and Y
                    Y = ThisPoint.Y
                    prePnt = ThisPart.getObject(point -1) # the vertex before this one
                    pstPnt = ThisPart.getObject(point +1) # the vertex after this one
                    # http://stackoverflow.com/questions/1211212/how-to-calculate-an-angle-from-three-points
                    # atan2(P2.y - P1.y, P2.x - P1.x) - atan2(P3.y - P1.y, P3.x - P1.x)
                    Angle = math.atan2(Y-prePnt.Y,X-prePnt.X) - math.atan2(pstPnt.Y - prePnt.Y,pstPnt.X - prePnt.X)
                    # I'm not sure if that's what you're after but relatively speaking that's the easy part
                    # to change if you're not happy with the angles. I am fairly sure the answer is in radians.
                    ICur.insertRow(((X,Y),Angle)) # write the point to the new shape file

I can't vouch for the maths on the last few lines but there's a fairly good skeleton there to work with.

The key is to get this vertex, the one before and the one after and through Pythagorean mathematics the interior angle can be determined.

Best Answer

Related Solutions

[GIS] Creating point at midpoint of polyline with ArcPy

[GIS] Calculating angle between polyline segments and assign this value to point

Related Question