Azimuth from last point to next point

I found an script on https://community.esri.com/blogs/dan_patterson/2016/09/01/distance-calculations-using-the-field-calculator but the script is giving the azimuth to the next Objectid point and what I really need is to get the azimuth from the last Point (based on point_m) placed on azimute1 Field and the azimuth to the next point placed on azimute2 using two codes on field calculator to do so. A personal Geodatabase with the point shape and the reference path line can be found here:Help.mdb.

``````x0 = 0.0
y0 = 0.0
angle = 0.0
def angle_between(shape):
global x0
global y0
x = shape.centroid.X
y = shape.centroid.Y
if x0 == 0.0 and y0 == 0.0:
x0 = x
y0 = y
return 0.0
radian = math.atan2((shape.centroid.Y - y0),(shape.centroid.X - x0))
x0 = x
y0 = y
return angle
angle_between(!Shape!)
``````

• Can you add a diagram of what you are trying to achieve to make the question clear? – Dan Mar 19 '18 at 0:14
• I´m trying to get the incoming azimuth and the outcoming azimuth from several points on a linear referenced feature – Luis Felipe Fragoso Mar 20 '18 at 1:12

This Python script should get you most of the way to what you're trying to achieve. It is a stand-alone script, not a Field Calculator expression. If you need to use it in ModelBuilder, you can add it to a custom tool and add it to ModelBuilder from there.

It takes the coordinates for the features and adds them into a Python list. This gives you the ability to access the properties of the points before and after the current one you are iterating on. Some basic coordinate geometry and trigonometry gives you the bearing and distance to each point. You may need to alter the output to give you exactly what you require. The script is currently set up to calculate the angles based on 0 degrees at north.

The output is directed to a table which you can then join to the input point feature class based on OBJECTID = LinkID.

``````import arcpy
import math

#---------------------------------------------------------------
def GetBearingDist(deltaX, deltaY):

if deltaX == 0.0 and deltaY == 0.0:
# points are concurrent
return [0.0, 0.0]
else:
# the angle between the points is one of the cardinal
# directions, no need for trigonometry.
if deltaX == 0.0 or deltaY == 0.0:
if deltaX == 0.0 and deltaY > 0.0:
angle = 0.0
elif deltaX == 0.0 and deltaY < 0.0:
angle = 180.0
elif deltaX > 0.0 and deltaY == 0.0:
angle = 90.0
elif deltaX < 0.0  and deltaY == 0.0:
angle = 270.0
else:
# calculate the angle in degrees and determine
# with 0 degrees being north.

if deltaX > 0.0 and deltaY > 0.0:
angle = 90.0 - abs(angle)
elif deltaX > 0.0 and deltaY < 0.0:
angle = 90.0 + abs(angle)
elif deltaX < 0.0 and deltaY < 0.0:
angle = 270.0 - abs(angle)
elif deltaX < 0.0 and deltaY > 0.0:
angle = 270.0 + abs(angle)

# calculate the distance while we are here
dist = math.sqrt(deltaX ** 2 + deltaY ** 2)

return [angle, dist]
#---------------------------------------------------------------

# Get input parameters from the user (hardcoded for testing purposes)
#aFC = arcpy.GetParameterAsText(0)
aFC = r"C:\temp\scratch.gdb\test_pnt"

#aOutTable = arcpy.GetParameterAsText(1)

# The output table must pre-exist with the following field definitions
# Back_Bearing -> Double
# Back_Distance -> Double
# Forward_Bearing -> Double
# Forward_Distance -> Double
# This could be automated, but is out of scope for this script.

# Build an array of the coordinates (in feature order)
# Add a SQL clause to the search cursor if you want the features in any other order
aCoordList = []

with arcpy.da.SearchCursor(aFC, ["OID@", "SHAPE@X", "SHAPE@Y"]) as aRows:
for aRow in aRows:
aCoordList.append([aRow[0], aRow[1], aRow[2]])

aMaxRows = len(aCoordList)

aOutput = []

for i in range(0, aMaxRows):
if i == 0:
# only calculate forward bearing for first record
dx = aCoordList[i + 1][1] - aCoordList[i][1]
dy = aCoordList[i + 1][2] - aCoordList[i][2]
brgdist = GetBearingDist(dx, dy)
aOutput.append([aCoordList[i][0], -999, -999, brgdist[0], brgdist[1]])
elif i == aMaxRows - 1:
# only calculate back bearing for last record
dx = aCoordList[i - 1][1] - aCoordList[i][1]
dy = aCoordList[i - 1][2] - aCoordList[i][2]
brgdist = GetBearingDist(dx, dy)
aOutput.append([aCoordList[i][0], brgdist[0], brgdist[1], -999, -999])
else:
# calculate forward and back bearings for all other records

# forward
dx = aCoordList[i + 1][1] - aCoordList[i][1]
dy = aCoordList[i + 1][2] - aCoordList[i][2]
brgdist_f = GetBearingDist(dx, dy)

# backward
dx = aCoordList[i - 1][1] - aCoordList[i][1]
dy = aCoordList[i - 1][2] - aCoordList[i][2]
brgdist_b = GetBearingDist(dx, dy)
aOutput.append([aCoordList[i][0], brgdist_b[0], brgdist_b[1], brgdist_f[0], brgdist_f[1]])

# write the output to a table so that it can be joined to the
# input feature class if required
with arcpy.da.InsertCursor(aOutTable, ["LinkID", "Back_Bearing", "Back_Distance", "Forward_Bearing", "Forward_Distance"]) as aIRows:
for k in aOutput:
aIRows.insertRow(k)

# finish
``````

Data from the output table for the diagram above.

``````[LinkID, Back_Bearing, Back_Distance, Forward_Bearing, Forward_Distance]
[1, -999, -999, 81.19320768041067, 2127.4117779304984]
[2, 261.19320768041064, 2127.4117779304984, 107.69088771442591, 1802.6438476383205]
[3, 287.6908877144259, 1802.6438476383205, 97.92692862922179, 1180.8898912440995]
[4, 277.9269286292218, 1180.8898912440995, 84.8055679201311, 1308.225016594544]
[5, 264.80556792013107, 1308.225016594544, 0.734523959818091, 1154.8959016067192]
[6, 180.7345239598181, 1154.8959016067192, 89.41537109097044, 1450.979334049593]
[7, 269.41537109097044, 1450.979334049593, 134.65691424063027, 1748.3235859907875]
[8, 314.6569142406303, 1748.3235859907875, 89.4327351833379, 1495.3925905719177]
[9, 269.4327351833379, 1495.3925905719177, 30.44357514873448, 1957.6966435050726]
[10, 210.44357514873448, 1957.6966435050726, 1.9091535067902043, 1777.6037365374682]
[11, 181.9091535067902, 1777.6037365374682, 270.0, 888.3085000002757]
[12, 90.0, 888.3085000002757, 206.56504831039555, 893.8432194864047]
[13, 26.565048310395547, 893.8432194864047, 229.51399125464334, 798.1055820384025]
[14, 49.513991254643344, 798.1055820384025, 269.01223580510293, 858.8259220405614]
[15, 89.01223580510293, 858.8259220405614, 324.6052066483075, 690.1483118337023]
[16, 144.60520664830753, 690.1483118337023, 7.594647593485973, 896.1697992178911]
[17, 187.59464759348597, 896.1697992178911, 46.50743669613549, 795.9054219035123]
[18, 226.5074366961355, 795.9054219035123, 71.32329730323522, 1109.5956812338072]
[19, 251.32329730323522, 1109.5956812338072, 58.626990496484865, 710.9551967988131]
[20, 238.62699049648486, 710.9551967988131, 29.001702431727985, 1557.3573559748613]
[21, 209.00170243172798, 1557.3573559748613, 21.889554140149585, 1787.0118803301573]
[22, 201.88955414014958, 1787.0118803301573, 90.0, 1347.2679000003263]
[23, 270.0, 1347.2679000003263, -999, -999]
``````

I thought I could answer your question with less code, but my answer is almost as long. It became just another way of doing it:

``````import arcpy
from itertools import tee

arcpy.env.workspace = r'C:\Arenden\Default.gdb' #Change to match your data
lines = 'Lines123'                              #Change to match your data
output_points = 'AnglePoints'                   #Change to match your data

def pairwise(iterable):
"s -> (s0,s1), (s1,s2), (s2, s3), ..."
a, b = tee(iterable)
next(b, None)
return zip(a, b)

def toangle_to_fromangle(e):
if e[0] is None:
return e
elif e[0]<=180:
e[0]+=180
else:
e[0]-=180
return e

#For each vertex list the angles
angles = []
with arcpy.da.SearchCursor(lines,'SHAPE@') as cursor:
for row in cursor:
for p in row[0]:
points = [arcpy.PointGeometry(a,arcpy.Describe(lines).spatialReference) for a in p]
for pair in pairwise(points):
angle = pair[0].angleAndDistanceTo(pair[1])[0]
if angle < 0:
angle += 360
angles.append(angle)

#Insert None for first and last vertex
angles.insert(0,None)
angles.insert(len(angles),None)

#Fore each vertex create pairs of from-to-angles
tofromangles = [list(p) for p in pairwise(angles)]
tofromangles = map(toangle_to_fromangle, tofromangles)

#Create a point feature class and insert geometries and from-to-angles
arcpy.CreateFeatureclass_management(out_path=arcpy.env.workspace, out_name=output_points,geometry_type='POINT',spatial_reference=arcpy.Describe(lines).spatialReference)
newfields = ['AzFrom','AzTo']
for f in newfields: