2

I have a polyline feature with many segments (polyline split at vertices).

I plan to record angle values to attribute table per each polyline segment and assign it to point features.

I need to calculate angle value between segments of a polyline (Get value calculated angle less 180 ). enter image description here.

Here almost is what I need: Calculating incidence angle of polyline?

5
  • 1
    Do you have any ArcObjects or arcpy ability? Mar 10, 2016 at 5:06
  • Lines and points.
    – Jannik
    Mar 10, 2016 at 10:54
  • What I'm asking is do you have any ability in scripting or programming? There might be a way to do this without scripting but it's very fiddly and interactive, it would be greatly preferable if you could do this in VB.net, C# or python, if so I can share some code that does approximately what you want to do. Mar 10, 2016 at 21:07
  • I think Python code more friendly for me. If you can, let me see method without scripting too.
    – Jannik
    Mar 11, 2016 at 1:48
  • To do this without scripting you must have an advanced license. Use Feature Vertices to points to generate the vertices as points then Generate Near Table which will record the angle (from baseline) between one and the other then you need to find the record prior and post for each point and calculate the difference in angle. To do this in scripting is much easier as you cursor through the line from second to second last and calculate the angle between previous and post, record as an XY table and generate points. If you can make some sense of python I'll try to throw something together. Mar 11, 2016 at 2:05

2 Answers 2

6

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.

2
  • +1 At first I thought you hadn't considered the first and last points but you did in qPointsRange. I haven't checked the math either, but it looks good to me.
    – Fezter
    Mar 11, 2016 at 3:19
  • 1
    @Michael Stimson I also was using math.atan function but when the true angle was about 175 degrees it could count it as 5 degrees. I turned to a bulky cosine function but it did rigth calculations. Apr 4, 2018 at 13:54
-1

Congrats to @Michael Stimson for the great script. I stumbled upon the thread when trying to tackle a similar problem and -after some test runs performed in a dataset of my own- I believe that the following corrections should be made:

# 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.degrees(math.atan2(Y-prePnt.Y,X-prePnt.X) - math.atan2(pstPnt.Y - prePnt.Y,pstPnt.X - prePnt.X))
                    Angle = math.degrees(math.atan2(prePnt.Y-Y,prePnt.X-X) - math.atan2(pstPnt.Y - Y,pstPnt.X - X))
                    if Angle > 0:
                        if Angle < 180:
                            Angle = Angle
                        else:
                            Angle = 360 -Angle
                    else:
                        if abs(Angle) < 180:
                            Angle = abs(Angle)
                        else:
                            Angle = 360 - abs(Angle)
                    ICur.insertRow(((X,Y),Angle)) # write the point to the new shape file

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.