7

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?

2
  • 1
    I checked azimuth between pt1 and pt2in field calculator: degrees( azimuth( make_point(30.30000000000000071, 60.60000000000000142), make_point(30.30129068361398836, 60.6006346117680792) ) ) It also results in 63.81727521235025
    – greenwd
    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
    – greenwd
    Feb 10, 2020 at 14:05

1 Answer 1

6

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.