Azimuth calculations using PyQGIS

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

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

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?

• I checked azimuth between `pt1` and `pt2`in field calculator: `degrees( azimuth( make_point(30.30000000000000071, 60.60000000000000142), make_point(30.30129068361398836, 60.6006346117680792) ) ) ` It also results in 63.81727521235025 Feb 5, 2020 at 14:04
• Also it is worth noticing that the length between pt1 and pt2 (`distance = da.measureLine(pt1, pt2)`) is 99.99999997504422. So the distance is measured properly, as opposed to the azimuth Feb 10, 2020 at 14:05

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

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
``````