Need some help with spherical trigonometry

[Trinopoty]

I need a formula to calculate the area of a triangle at radius r given 3 points (λ1, Φ1), (λ2, Φ2), (λ3, Φ3).

I have looked at https://en.wikipedia.org/wiki/Spherical_trigonometry and came up with an implementation in python but the sum of the angles always came up to slightly less than 2π making E negative.

I don't know if this is a problem with floating point precision or some other issue.

For the record, my points are very close to each other and are latitude and longitude and suppose r to be the mean radius of earth at sea level.

 

[Franck]

You know calculating spherical trigonometry is not the way to calculate earth surface. You will be off by (if i recall correctly) +12% to -5%

[Trinopoty]

3 hours ago, Franck said:

You know calculating spherical trigonometry is not the way to calculate earth surface. You will be off by (if i recall correctly) +12% to -5%

Well. I heard (read?) that it's the most accurate. I have 3 different result for the same coordinates and I'm trying to verify which one is the correct one using spherical trigonometry method.

[0x21]

Use Girard's theorem.

 

In the case that your angles sum to ≈ 2π, then E should be positive and ≈ π.

[Trinopoty]

4 hours ago, Serjeant said:

Use Girard's theorem.

 

In the case that your angles sum to ≈ 2π, then E should be positive and ≈ π.

I am using Girard's theorem and the sum of angles are < π which should not be possible. My primary suspect is floating point error.

[wanderingfool2]

On 3/20/2020 at 1:06 AM, Trinopoty said:

For the record, my points are very close to each other and are latitude and longitude and suppose r to be the mean radius of earth at sea level

How close is very close?  (Do you have the numbers you are using/examples of how close, or the code you use).

 

If the numbers are really close, then yes it might be floating point rounding...python gives 15-17 (base 10) digits of accuracy so keep that in mind...if you positions are close enough that only the first 13 digits match, for example, then you will get quite a bad result

[0x21]

You could use L’Huilier’s Theorem and do everything in terms of side lengths, which are trivial to calculate with polar coordinates.

[straight_stewie]

On 3/20/2020 at 3:06 AM, Trinopoty said:

but the sum of the angles always came up to slightly less than 2π making E negative.

Is it guaranteed that the points define a triangle and not an open polygon?

By definition, if the sum of the three angles is less than 180, you have an open polygon.
 

21 hours ago, Trinopoty said:

My primary suspect is floating point error.

Have you worked the failure case(s) by hand to make sure that they have valid solutions?

[Trinopoty]

18 hours ago, wanderingfool2 said:

How close is very close?

Area of triangle as calculated by QGIS is about 600 m sq.

 

3 hours ago, straight_stewie said:

Is it guaranteed that the points define a triangle and not an open polygon?

Yes the points form a triangle not a polygon.

 

3 hours ago, straight_stewie said:

Have you worked the failure case(s) by hand to make sure that they have valid solutions?

QGIS is able to calculate valid values. So are other libraries such as openlayers. But the values vary wildly and that's what I'm trying to check.

 

17 hours ago, Serjeant said:

You could use L’Huilier’s Theorem

I'll try that.

[Trinopoty]

My code:

import math
from math import sin
from math import cos
from math import asin
from math import acos
from math import sqrt

import json


def sinsq(θ):
    return ((1 - cos(2 * θ)) / 2)


def _Δσ1(p1, p2):
    λ1 = p1[0]
    Φ1 = p1[1]
    λ2 = p2[0]
    Φ2 = p2[1]
    
    return acos((sin(Φ1) * sin(Φ2)) + (cos(Φ1) * cos(Φ2) * cos(λ1 - λ2)))


def _Δσ2(p1, p2):
    λ1 = p1[0]
    Φ1 = p1[1]
    λ2 = p2[0]
    Φ2 = p2[1]
    
    ΔΦ = abs(Φ1 - Φ2)
    Δλ = abs(λ1 - λ2)
    
    return 2 * asin(sqrt(sinsq(ΔΦ / 2) + (cos(Φ1) * cos(Φ2) * sinsq(Δλ / 2))))


def triangle_area1(r, pA, pB, pC):
    a1 = _Δσ1(pB, pC)
    b1 = _Δσ1(pA, pC)
    c1 = _Δσ1(pA, pB)
    
    α1 = acos((cos(a1) - (cos(b1) * cos(c1))) / (sin(b1) * sin(c1)))
    β1 = acos((cos(b1) - (cos(a1) * cos(c1))) / (sin(a1) * sin(c1)))
    Ɣ1 = acos((cos(c1) - (cos(a1) * cos(b1))) / (sin(a1) * sin(b1)))
    e1 = α1 + β1 + Ɣ1 - math.pi
    
    return r * r * e1


def triangle_area2(r, pA, pB, pC):
    a2 = _Δσ2(pB, pC)
    b2 = _Δσ2(pA, pC)
    c2 = _Δσ2(pA, pB)
    
    α2 = acos((cos(a2) - (cos(b2) * cos(c2))) / (sin(b2) * sin(c2)))
    β2 = acos((cos(b2) - (cos(a2) * cos(c2))) / (sin(a2) * sin(c2)))
    Ɣ2 = acos((cos(c2) - (cos(a2) * cos(b2))) / (sin(a2) * sin(b2)))
    e2 = α2 + β2 + Ɣ2 - math.pi
    
    return r * r * e2


def main():
    # Area by QGIS: 72668327810.760 m²
    triangles = [[[75.6298828125,21.22794190505815],[77.16796875,18.729501999072138],[79.991455078125,22.339914425562032]]]
    triangles = [[[math.radians(v) for v in p] for p in triangle] for triangle in triangles]
    
    total_area = 0
    for triangle in triangles:
        area = triangle_area1(6378137, triangle[0], triangle[1], triangle[2])
        total_area = total_area + area
    print('Total: ', total_area)


if __name__ == '__main__':
    main()