Azimuth calculations using PyQGIS

Question

According to the QgsDistanceArea docs, an azimuth for computeSpheroidProject should be:

in radians, clockwise from North

So, to build a new point with computeSpheroidProject I use math.radians(azimuth)

pt1 = QgsPointXY(30.30, 60.60)
distance = 100
azimuth = 45

da = QgsDistanceArea()
da.setEllipsoid('WGS84')
pt2 = da.computeSpheroidProject(pt1, distance, math.radians(azimuth))

print(pt1, pt2)
>>><QgsPointXY: POINT(30.30000000000000071 60.60000000000000142)> 
<QgsPointXY: POINT(30.30129068361398836 60.6006346117680792)>

But when needed to find an azimuth between pt1 and pt2 it doesn't return the original value (45):

azimuth = pt1.azimuth(pt2)
print(azimuth)
>>>63.81727521235025

Looks like I haven't quite figured out azimuth calculations. Why results are different?

Answer 1

The results are different because QgsPointXY::azimuth does not consider the ellipsoid. It instead calculates the azimuth on a cartesian plane by using the formula or in python terms math.atan2(pt2.y() - pt1.y(), pt2.x() - pt1.x()). I'm not even sure if you can call the result of such computation an azimuth as the method suggests because in general azimuth seems to be an angular measurement in a spherical coordinate system not on a plane.

computeSpheroidProject however does consider that the earth is not a flat plane and performs a geodesic ("Great Circle") computation of the azimuth. For the backwards computation you should therefore also calculate with respect to the ellipsoid by using QgsDistanceArea again. The function for that is QgsDistanceArea::bearing which returns a result in radians.

pt1 = QgsPointXY(30.30, 60.60)
distance = 100
azimuth = 45

da = QgsDistanceArea()
da.setEllipsoid('WGS84')
pt2 = da.computeSpheroidProject(pt1, distance, math.radians(azimuth))

azimuth = math.degrees(da.bearing(pt1, pt2))
print(azimuth) # 44.99999999951381

or for computing in cartesian coordinates use the methods from QgsPointXY class for each direction of the computation:

pt2 = pt1.project(distance, azimuth)
azimuth = pt1.azimuth(pt2)
print(azimuth) # 45.0