Calculating azimuth between two points using their projected coordinates with Python

Question

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?

Answer 1

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).

Answer 2

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
network = gpd.read_file(r'path/to/file.gpkg')
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

• Andrew Tice: Things are usually a threat when trig is involved ;)
• i.i.k.: @AndrewTice Hehe, I adjusted it ;) Thanks for the hint