4

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?

2 Answers 2

5

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

0
3

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
2
  • Things are usually a threat when trig is involved ;) Apr 23 at 22:13
  • 1
    @AndrewTice Hehe, I adjusted it ;) Thanks for the hint
    – i.i.k.
    Apr 24 at 7:21

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.