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:
The result doesn't look like the two are at perpendicular angles:
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.
Your computation is correct. However, a 288.7° angle with respect to the X axis should look like this:
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:
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°