Search code examples
geopandasshapefileqgisshapelyazimuth

calculate linestring endpoints from point, perpendicular to a given angle

I am trying to use shapely and trigonometry to calculate endpoints of a line that passes through a point, perpendicularly to a given angle. For example, if my angle is 90, the line should have an azimuth of 0 or 180 degrees. The endpoints order doesn't really matter so both values would work fine.

Here's a minimum reproducible example of the code:

import math
from shapely import Point, LineString
import geopandas as gpd


def calculate_endpoints(start_point, azimuth_deg, distance = 1000):

    x, y = x, y = start_point.coords[0]
    azimuth_rad = math.radians(azimuth_deg)
    
    right_azimuth_rad = (azimuth_rad + math.pi / 2) % (2 * math.pi)
    left_azimuth_rad = (azimuth_rad - math.pi / 2) % (2 * math.pi)
    
    dx_right = distance * math.cos(right_azimuth_rad)
    dy_right = distance * math.sin(right_azimuth_rad)
    right_endpoint = (x + dx_right, y + dy_right)
    
    dx_left = distance * math.cos(left_azimuth_rad)
    dy_left = distance * math.sin(left_azimuth_rad)
    left_endpoint = (x + dx_left, y + dy_left)
    
    return LineString([right_endpoint, left_endpoint])

start_point = Point(8424, 1385)
azimuth = 288.7

If I save the results to a shapefile:

gpd.GeoDataFrame( geometry=[calculate_endpoints(start_point, azimuth, 1000)]).to_file("line.shp")

gpd.GeoDataFrame( geometry=[start_point]).to_file("point.shp")

But then in QGIS:

  • Import the point shapefile and change the point symbology to line, and rotate it by 288.7 degrees enter image description here
  • Add the line shapefile

The result doesn't look like the two are at perpendicular angles:

enter image description here

But if I calculate the line's azimuth,

x1 = calculate_endpoints(start_point, azimuth, 1000).coords.xy[0][0]
y1 = calculate_endpoints(start_point, azimuth, 1000).coords.xy[0][1]
x2 = calculate_endpoints(start_point, azimuth, 1000).coords.xy[1][0]
y2 = calculate_endpoints(start_point, azimuth, 1000).coords.xy[1][1]

dx = x2-x1
dy = y2-y1
math.degrees(math.atan2(dy, dx))

The result is -140, which is roughly perpendicular to the original angle? So am I visualizing this wrong in QGIS, or is there something wrong with the calculations? Any help is appreciated.

Solution

  • Your computation is correct. However, a 288.7° angle with respect to the X axis should look like this:

    enter image description here

    According to the QGIS documentation, the marker for a point can be a vertical line. In addition, rotations are clockwise with respect to North in QGIS, see this post so the result is this:

    enter image description here

    If you want the marker to be correct, you need to change your marker rotation angle in QGIS to 90 + (360-277.8) = 161.3°