A few days ago I wrote about how you could calculate the area of a spherical triangle by measuring the angles at its vertices. As was pointed out in the comments, there’s a much more direct way to find the area.

Folke Eriksson gives the following formula for area in [1]. If **a**, **b**, and **c** are three unit vectors, the spherical triangle formed by their vertices has area *E* where

Let’s try the formula out with an example and compare the area to that of a flat triangle. We’ll pick three cities as our three points, assuming the earth is a unit sphere.

Let **a** = Austin, Texas, **b** = Boston, Massachusetts, and **c** = Calgary, Alberta. The spherical coordinates of each point will be (1, θ, φ) where θ is longitude and φ is π/2 minus latitude.

from numpy import * def latlong_to_cartesian(latlong): phi = deg2rad( 90 - latlong[0] ) theta = deg2rad( latlong[1] ) x = sin(theta) * cos(phi) y = sin(theta) * sin(phi) z = cos(theta) return array([x,y,z]) def spherical_area(a, b, c): t = abs( inner(a, cross(b, c) ) ) t /= 1 + inner(a,b) + inner(b,c) + inner(a,c) return 2*arctan(t) def flat_area(a, b, c): x = cross(a - b, c - b) return 0.5*linalg.norm(x) # Latitude and longitude austin = (30.2672, 97.7431) boston = (42.3584, 71.0598) calgary = (51.0501, 114.0853) a = latlong_to_cartesian(austin) b = latlong_to_cartesian(boston) c = latlong_to_cartesian(calgary) round = spherical_area(a, b, c) flat = flat_area(a, b, c) print(round, flat, round/flat)

This shows that our spherical triangle has area 0.1154 and the corresponding Euclidean triangle slicing through the earth would have area 0.1093, the former being 5.57% bigger.

Let’s repeat our exercise with a smaller triangle. We’ll replace Boston and Calgary with Texas cities Boerne and Carrollton.

boerne = (29.7947, 98.7320) carrollton = (32.9746, 96.8899)

Now we get area 0.000345 in both cases. The triangle is so small relative to the size of the earth that the effect of the earth’s curvature is negligible.

In fact the ratio of the two computed areas is 0.9999, i.e. the spherical triangle has slightly less area. The triangle is so flat that numerical error has a bigger effect than the curvature of the earth.

How could we change our code to get more accurate results for small triangles? That would be an interesting problem. Maybe I’ll look into that in another post.

[1] “On the measure of solid angles”, F. Eriksson, Math. Mag. 63 (1990) 184–187.

Use more precise coordinates. Such as state plane.

More precise data would give more accurate results, taking into account deviations of the earth from being perfectly spherical.

But that wouldn’t address the numerical accuracy problem.

I think the phi calculation in latlong_to_cartesian might be wrong. Pretty sure the pi/2 subtraction should be outside of the deg2rad call.

Thanks, Guillaume. You’re right. Interestingly, it doesn’t make much numerical difference in this case.