Azimuth calculation in QGIS

Question

In QGIS 3.10 I am calculating the azimuth between two points whose position is defined in WGS84 coordinates. If I Open Attribute Table -> Open Field Calculator and ask for

degrees(
    azimuth(
        make_point(6.020323,47.248690),
        make_point(6.083703,47.242157)
        )
    )

I get 95.885° degrees which is inconsistent with the manual calculation:

atan2(sin(delta_lon)*cos(lat2), cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(delta_lon))

or with the result provided by the French Geographic Institute (Géoportail) website.

However, if following Azimuth calculation, extending lines I run in the Python Console:

p1 = QgsPointXY(6.020323,47.248690)
p2 = QgsPointXY(6.083703,47.242157)
d = QgsDistanceArea()
d.setEllipsoid('WGS84')
d.bearing(p1,p2)*180/3.1415926535

then I get the correct result of 98.584. Could someone enlighten me as to why the first method provides an erroneous result?

Answer 1

The azimuth() expression function calculates the Cartesian angle y-axis-based clockwise in the Cartesian plane of the coordinate reference system in which are expressed the coordinates of the given points.

To obtain a closer result (True-North-based clockwise, loxodromic) using a Cartesian angle calculation, you shouldn't use a geographical coordinate reference system, but e.g. a conformal projected one, as a proper UTM projection, and then correct the result taking into account the grid convergence at the first point (see Calculating convergence angle for point data in QGIS).

The QgsDistanceArea.bearing() PyQGIS function, instead, calculates the geodesic bearing angle in radians True-North-based clockwise on the specified ellipsoid (orthodromic) from the given WGS 84 geographical coordinates, or from the coordinates expressed in another coordinate reference system just setting it with the QgsDistanceArea function setSourceCrs()

UPDATE: since QGIS 3.34, it is also available the bearing() expression function, which is based on the above mentioned QgsDistanceArea.bearing() function and returns the north-based bearing as the angle in radians measured clockwise on the ellipsoid from the vertical on point_a to point_b in the specified CRS and on the specified ellipsoid.