# Calculating azimuth between two points using their projected coordinates with Python

First of all, the goal is to create small lines from points at each end of a LineString. However, for it a step is to calculate the azimuth for the first and second, and second last and last vertex.

First I import a GeoPackage (CRS is EPSG:25833) as GeoPandas `gdf` and extract a single line.

``````network = gpd.read_file(r'path/to/file.gpkg')
singleLineFromNetwork = network.iloc[558].geometry
singleLineFromNetwork
>>> shapely.geometry.linestring.LineString
``````

I then create point geometries based on the coordinates of the first, second, second last and last vertices.

``````first = shapely.geometry.Point(singleLineFromNetwork.coords[0])
second = shapely.geometry.Point(singleLineFromNetwork.coords[1])

last = shapely.geometry.Point(singleLineFromNetwork.coords[len(line.coords.xy[0])-1])
secondlast = shapely.geometry.Point(singleLineFromNetwork.coords[len(line.coords.xy[0])-2])
``````

How can I now calculate the azimuth between the first and second geometry and from the last and second last points?

Is there some function in shapely or GeoPandas, or how can I do it with e.g. math libraries?

Since your coordinates are on a Cartesian plane you can use the trigonometry that you learned in high school. Your coordinates are a right angled triangle with sides of `dx=Math.abs(x1-x2)` and `dy=Math.abs(y1-y2)` and since you don't necessarily know the hypotenuse you can use `tan(theta) = dx/dy` to give you the angle of the line to the Y-axis (you may need to add pi/2 to allow for which direction you are going in).

Based on this thread (How is it possible to calculate azimuth from a single lat lon cordinate with x, y offsets in meters?) I was able to calculate the azimuth using the math library in Python (thanks @Hossein). The script is as follows :

``````def calculate_azimuth(ax, ay, bx, by):
"""Computes the bearing in degrees from the point A(a1,a2) to
the point B(b1,b2). Note that A and B are given in terms of
screen coordinates.

Args:
ax (int): The x-coordinate of the first point defining a line.
ay (int): The y-coordinate of the first point defining a line.
bx (int): The x-coordinate of the second point defining a line.
by (int): The y-coordinate of the second point defining a line.

Returns:
float: bearing in degrees
"""

TWO_PI = math.pi * 2
# if (a1 = b1 and a2 = b2) throw an error
theta = math.atan2(bx - ax, ay - by)
if (theta < 0.0):
theta += TWO_PI
return math.degrees(theta)

``````

Integrating this with my problem above:

``````import math, shapely
import geopandas as gpd

singleLineFromNetwork = network.iloc[558].geometry

first = shapely.geometry.Point(singleLineFromNetwork.coords[0])
second = shapely.geometry.Point(singleLineFromNetwork.coords[1])

last = shapely.geometry.Point(singleLineFromNetwork.coords[len(line.coords.xy[0])-1])
secondlast = shapely.geometry.Point(singleLineFromNetwork.coords[len(line.coords.xy[0])-2])

azimuthStartOfLine = calculate_azimuth(first[0], first[1], second[0], second[1])
azimuthEndOfLine = calculate_azimuth(secondlast[0], secondlast[1], last[0], last[1])
``````

For example:

``````calculate_azimuth(415162.39785234944, 5928033.199373602, 415172.0, 5928036.0)

>>> 106.26020470823194
``````
• Things are usually a threat when trig is involved ;) Apr 23 at 22:13
• @AndrewTice Hehe, I adjusted it ;) Thanks for the hint Apr 24 at 7:21