# Spherical trig, Research Triangle, and Mathematica

This post will look at the triangle behind North Carolina’s Research Triangle using Mathematica’s geographic functions.

## Spherical triangles

A spherical triangle is a triangle drawn on the surface of a sphere. It has three vertices, given by points on the sphere, and three sides. The sides of the triangle are portions of great circles running between two vertices. A great circle is a circle of maximum radius, a circle with the same center as the sphere.

An interesting aspect of spherical geometry is that both the sides and angles of a spherical triangle are angles. Because the sides of a spherical triangle are arcs, they have angular measure, the angle formed by connecting each vertex to the center of the sphere. The arc length of a side is its angular measure times the radius of the sphere.

Denote the three vertices by AB, and C. Denote the side opposite A by a, etc. Denote the angles at AB, and C by α, β, and γ respectively.

## Research Triangle

Research Triangle is a (spherical!) triangle with vertices formed by Duke University, North Carolina State University, and University of North Carolina at Chapel Hill.

(That was the origin of the name, though it’s now applied more loosely to the general area around these three universities.)

We’ll take as our vertices

• A = UNC Chapel Hill (35.9046 N, 79.0468 W)
• B = Duke University in Durham (36.0011 N, 78.9389 W),
• C = NCSU in Raleigh (35.7872 N, 78.6705 W) ## Mathematica

We’ll illustrate several features of Mathematica using the spherical triangle corresponding to Research Triangle.

### Map

The map above was produced with the following Mathematica code.

    ptA = GeoPosition[{35.9046, -79.0468}]
ptB = GeoPosition[{36.0011, -78.9389}]
ptC = GeoPosition[{35.7872, -78.6705}]

GeoGraphics[{Red, PointSize[Large],
Point[ptA], Point[ptB], Point[ptC]},
GeoScaleBar -> "Imperial",
GeoRange -> 29000]


Note that longitude is in degrees east, so the longitudes above are negative.

### Distance

To find the distance between two locations, you can use the function GeoDistance. For example,

    GeoDistance[ptA, ptB]


tells me that the distance between UNC and Duke is 8.99185 miles. I assume it displays miles by default based on my location in the US, though the GeoRange option above defaults to meters. You can make the system of units explicit. For example

    GeoDistance[ptA, ptB, UnitSystem -> "Metric"]


returns 14.741 km.

If we want to find the length of the side c in degrees, we can use the Earth’s radius.

    r = PlanetData["Earth", "Radius"]
c = GeoDistance[ptA, ptB] 360 / (2 Pi r)


This shows that c is 0.13014°. Calculating the other sides similarly shows a = 0.30504° and b = 0.32739°.

### Angles

Calling GeoDirection[ptA, ptB] returns 42.2432°, which says we need to head at an angle of about 42° to walk from UNC to Duke.

The code

    GeoDirection[ptA, ptB] - GeoDirection[ptA, ptC]

shows that the angle α is 68.6128°. (The code returns the negative of this angle because the angle is clockwise.) Similarly we find β = 87.9808° and γ = 23.4068.

The angles add up to 180° only because our triangle is small compared to the earth’s total area. The actual sum should be slightly more than 180°, but we’ve not retained enough precision to detect the difference. In general the “spherical excess,” i.e. the amount by which the sum of the angles exceed 180°, is proportional to the area of the triangle.

# Ellipsoid distance on Earth To first approximation, Earth is a sphere. But it bulges at the equator, and to second approximation, Earth is an oblate spheroid. Earth is not exactly an oblate spheroid either, but the error in the oblate spheroid model is about 100x smaller than the error in the spherical model.

Finding the distance between two points on a sphere is fairly simple. Here’s a calculator to compute the distance, and here’s a derivation of the formula used in the calculator.

Finding the distance between two points on an ellipsoid is much more complicated. (A spheroid is a kind of ellipsoid.) Wikipedia gives a description of Vincenty’s algorithm for finding the distance between two points on Earth using an oblate spheroid model (specifically WGS-84). I’ll include a Python implementation below.

## Comparison with spherical distance

How much difference does it make when you calculate difference on an oblate spheroid rather than a sphere? To address that question I looked at the coordinates of several cities around the world using the CityData function in Mathematica. Latitude is in degrees north of the equator and longitude is in degrees east of the prime meridian.

    |--------------+--------+---------|
| City         |    Lat |    Long |
|--------------+--------+---------|
| Houston      |  29.78 |  -95.39 |
| Caracas      |  10.54 |  -66.93 |
| London       |  51.50 |   -0.12 |
| Tokyo        |  35.67 |  139.77 |
| Delhi        |  28.67 |   77.21 |
| Honolulu     |  21.31 | -157.83 |
| Sao Paulo    | -23.53 |  -46.63 |
| New York     |  40.66 |  -73.94 |
| Los Angeles  |  34.02 | -118.41 |
| Cape Town    | -33.93 |   18.46 |
| Sydney       | -33.87 |  151.21 |
| Tromsø       |  69.66 |   18.94 |
| Singapore    |   1.30 |  103.85 |
|--------------+--------+---------|


Here are the error extremes.

The spherical model underestimates the distance from London to Tokyo by 12.88 km, and it overestimates the distance from London to Cape Town by 45.40 km.

The relative error is most negative for London to New York (-0.157%) and most positive for Tokyo to Sidney (0.545%).

Update: The comparison above used the same radius for both the spherical and ellipsoidal models. As suggested in the comments, a better comparison would use the mean radius (average of equatorial and polar radii) in the spherical model rather than the equatorial radius.

With that change the most negative absolute error is between Tokyo and New York at -30,038 m. The most negative relative error is between London and New York at -0.324%.

The largest positive error, absolute and relative, is between Tokyo and Sydney. The absolute error is 29,289 m and the relative error is 0.376%.

## Python implementation

The code below is a direct implementation of the equations in the Wikipedia article.

Note that longitude and latitude below are assumed to be in radians. You can convert from degrees to radians with SciPy’s deg2rad function.

from scipy import sin, cos, tan, arctan, arctan2, arccos, pi

def spherical_distance(lat1, long1, lat2, long2):
phi1 = 0.5*pi - lat1
phi2 = 0.5*pi - lat2
r = 0.5*(6378137 + 6356752) # mean radius in meters
t = sin(phi1)*sin(phi2)*cos(long1-long2) + cos(phi1)*cos(phi2)
return r * arccos(t)

def ellipsoidal_distance(lat1, long1, lat2, long2):

a = 6378137.0 # equatorial radius in meters
f = 1/298.257223563 # ellipsoid flattening
b = (1 - f)*a
tolerance = 1e-11 # to stop iteration

phi1, phi2 = lat1, lat2
U1 = arctan((1-f)*tan(phi1))
U2 = arctan((1-f)*tan(phi2))
L1, L2 = long1, long2
L = L2 - L1

lambda_old = L + 0

while True:

t = (cos(U2)*sin(lambda_old))**2
t += (cos(U1)*sin(U2) - sin(U1)*cos(U2)*cos(lambda_old))**2
sin_sigma = t**0.5
cos_sigma = sin(U1)*sin(U2) + cos(U1)*cos(U2)*cos(lambda_old)
sigma = arctan2(sin_sigma, cos_sigma)

sin_alpha = cos(U1)*cos(U2)*sin(lambda_old) / sin_sigma
cos_sq_alpha = 1 - sin_alpha**2
cos_2sigma_m = cos_sigma - 2*sin(U1)*sin(U2)/cos_sq_alpha
C = f*cos_sq_alpha*(4 + f*(4-3*cos_sq_alpha))/16

t = sigma + C*sin_sigma*(cos_2sigma_m + C*cos_sigma*(-1 + 2*cos_2sigma_m**2))
lambda_new = L + (1 - C)*f*sin_alpha*t
if abs(lambda_new - lambda_old) <= tolerance:
break
else:
lambda_old = lambda_new

u2 = cos_sq_alpha*((a**2 - b**2)/b**2)
A = 1 + (u2/16384)*(4096 + u2*(-768+u2*(320 - 175*u2)))
B = (u2/1024)*(256 + u2*(-128 + u2*(74 - 47*u2)))
t = cos_2sigma_m + 0.25*B*(cos_sigma*(-1 + 2*cos_2sigma_m**2))
t -= (B/6)*cos_2sigma_m*(-3 + 4*sin_sigma**2)*(-3 + 4*cos_2sigma_m**2)
delta_sigma = B * sin_sigma * t
s = b*A*(sigma - delta_sigma)

return s


# Projecting the globe onto regular solids

I was playing around with some geographic features of Mathematica this morning and ran across an interesting example in the documentation for the PolyhedronProjection function given here. Here’s what you get when you project a map of the earth onto each of the five regular (Platonic) solids.     ## How the images were made

At first I right-clicked on each image and saved as graphic files. This produced low quality images, even when I saved as SVG. I got better resolution by using the Export command and specifying the ImageSize and ImageResolution options.

The default view of the tetrahedron was face-on, so it looked like a flat triangle. By changing the ViewPoint I was able to get something that’s more clearly three dimensional.

By the way, you can use PolyhedronProjection to project onto more exotic polyhedra than the Platonic solids. For example,

    Export["rhomb.png",
PolyhedronProjection["ParagyrateDiminishedRhombicosidodecahedron"],
ImageResolution -> 72,
ImageSize -> Large]


produces this: # Equal Earth map projection

There’s no perfectly satisfying way to map the globe on to a flat surface. Every projection has its advantages and disadvantages. The Mercator projection, for example, is much maligned for the way it distorts area, but it has the property that lines of constant bearing correspond to straight lines on the map. Obviously this is convenient if you’re sailing without GPS. But for contemporary use, say in a classroom, minimizing area distortion is often a higher priority than keeping bearing lines straight.

Bojan Šavrič, Tom Patterson, and Bernhard Jenny have developed a new map projection called Equal Earth that nicely balances several competing criteria, including aesthetics. The Equal Earth projection satisfies the following mathematical criteria:

1. Equal area: The area of a region on the globe is proportional to the area of its projection.
2. Straight parallels: Lines of equal latitude project on to horizontal lines.
3. Regularly distributed meridians: At a given latitude, the spacing between lines of longitude are equal.
4. Bilateral symmetry: The projection is symmetric with respect to the x-axis (equator) and y-axis (prime meridian).

Here are the equations for the Equal Earth projection. Here λ and φ are longitude and latitude respectively. The parametric latitude θ is introduced for convenience and is simply a rescaling of latitude φ.

The parameters ai are given below.

    a_1 =  1.340264
a_2 = -0.081106
a_3 =  0.000893
a_4 =  0.003796


You can find more in their paper: Bojan Šavrič, Tom Patterson, and Bernhard Jenny. The Equal Earth map projection. International Journal of Geographical Information Science. https://doi.org/10.1080/13658816.2018.1504949.

Related posts

# Misplacing a continent

There are many conventions for describing points on a sphere. For example, does latitude zero refer to the North Pole or the equator? Mathematicians tend to prefer the former and geoscientists the latter. There are also varying conventions for longitude.

Volker Michel describes this clash of conventions colorfully in his book on constructive approximation.

Many mathematicians have faced weird jigsaw puzzles with misplaced continents after using a data set from a geoscientist. If you ever get such figures, too, or if you are, for example, desperately searching South America in a data set but cannot find it, remember the remark you have just read to solve your problem.

Related posts:

# Ellipsoid surface area

How much difference does the earth’s equatorial bulge make in its surface area?

To first approximation, the earth is a sphere. The next step in sophistication is to model the earth as an ellipsoid.

The surface area of an ellipsoid with semi-axes abc is where and The functions E and F are incomplete elliptic integrals and implemented in SciPy as ellipeinc and ellipkinc. Note that the SciPy functions take m as their second argument rather its square root k.

For the earth, a = b and so m = 1.

The following Python code computes the ratio of earth’s surface area as an ellipsoid to its area as a sphere.

from scipy import pi, sin, cos, arccos
from scipy.special import ellipkinc, ellipeinc

# values in meters based on GRS 80
# http://en.wikipedia.org/wiki/GRS_80

a = b = equatorial_radius

phi = arccos(c/a)
# in general, m = (a**2 * (b**2 - c**2)) / (b**2 * (a**2 - c**2))
m = 1

temp = ellipeinc(phi, m)*sin(phi)**2 + ellipkinc(phi, m)*cos(phi)**2
ellipsoid_area = 2*pi*(c**2 + a*b*temp/sin(phi))

# sphere with radius equal to average of polar and equatorial
r = 0.5*(a+c)
sphere_area = 4*pi*r**2

print(ellipsoid_area/sphere_area)


This shows that the ellipsoid model leads to 0.112% more surface area relative to a sphere.

Source: See equation 19.33.2 here.

Update: It was suggested in the comments that it would be better to compare the ellipsoid area to that of a sphere of the same volume. So instead of using the average of the polar and equatorial radii, one would take the geometric mean of the polar radius and two copies of the equatorial radius. Using that radius, the ellipsoid has 0.0002% more area than the sphere.

For daily posts on analysis, follow @AnalysisFact on Twitter. # A book I’d like (someone) to write

Here’s an idea I had for a book. Maybe someone has already written it. If you know of such a book, please let me know.

Differential geometry has a huge ratio of definitions to theorems. It seems like you do nothing but study definitions for a semester or two in preparation for proving something later. It’s easy to lose sight of the geometry. I’d like to see a book that is a concrete complement more typical abstract books.

My suggestion is for someone to write a book that goes through a standard differential geometry book, like Spivak’s, and compute everything for a small number of example manifolds: at least a sphere and an ellipsoid, maybe a torus. The book would first go through everything on a sphere where things are simplest, then generalize to an ellipsoid. There would be a lot of applications to geodesy: to first approximation the earth is a sphere, to second approximation it is an ellipsoid.

Sometimes a calculation, such as arc length, is very simple on a sphere. It can be done just using trig. Then the analogous calculation on an ellipsoid is much harder. It is complicated enough to illustrate the machinery of differential geometry. However, we know the answers shouldn’t be much different from those for a sphere, so we have a way to see whether the results are reasonable. This book would not shy away from computational difficulties.

I imagine this book would have lots of illustrations. It might even come with physical models, such as a globe with an exaggerated equatorial bulge. The idea is to be as tangible as other books are abstract.

I don’t plan to write this book, at least not any time soon. Maybe if my consulting goes well I would have the time to work on it in the future, but now is not the time for me to write a book. In the mean time, if someone wants to scoop my idea, please do!

* * *

Here are a couple other book ideas I’ve blogged about: R: The Good Parts and a rigorous elementary statistics text.

And here are some posts on geodesy: What is the shape of the earth? and Latitude doesn’t exactly mean what I thought.

And finally a few posts on spherical geometry: Napier’s mnemonic, The Sydney Opera House, and Mercator projection.  # Spotting sensitivity in an equation

The new book Heavenly Mathematics describes in the first chapter how the medieval scholar Abū Rayḥān al-Bīrūnī calculated the earth’s radius. The derivation itself is interesting, but here I want to expand on a parenthetical remark about the calculation.

The earth’s radius r can be found by solving the following equation. The constant in the denominator comes from a mountain which is 305.1 meters tall. The angle θ is known to be 34 minutes, i.e. 34/60 degrees. Here is the remark that caught my eye as someone more interested in numerical analysis than trigonometry:

There is a delicate matter hidden in this solution however: a minute change in the value of θ results in a large change in the value of r.

How can you tell that the solution is sensitive to changes (i.e. measurement errors) in θ? That doesn’t seem obvious.

Think of r as a function of θ and differentiate both sides of the equation with respect to θ. We’ll convert θ to radians because that’s what we do. (Explanation at the bottom of this post.) We get or Now let’s get a feel for the size of the terms in this equation. θ is approximately 0.01 radians, and so sin θ is approximately 0.01 as well. (See explanation here.) The radius of the earth is about 6.4 million meters. So the right side of the equation above is about 1.3 billion meters, i.e. it’s big.

A tiny increase in θ leads to a large decrease in r. For example, if our measurement of θ increased by 1%, from 0.01 to 0.0101, our measurement of the earth’s radius would decrease by 130,000 meters.

I’d like to point out a couple things about this analysis. First, it shows how it can be useful to think of constants as variables. After measuring θ we could think that we know its value with certainty and treat it as a constant. But a more sophisticated analysis takes into account that while θ might not change, our measurement of θ has changed from the true value.

Second, we used the radius of the earth to determine how sensitive our estimate of the earth’s radius is to changes in θ. Isn’t that circular reasoning? Not really. We can use a very crude estimate of the earth’s radius to estimate how sensitive a new estimate is to changes in its parameters. You always have some idea how big a value is before you measure it. If you want to measure the distance to the moon, you know not to pick up a yard stick. # Approximating Earth as a sphere

Isaac Newton suggested in 1687 that the earth is not a perfectly round sphere but rather an ellipsoid, and he was right. But since our planet is roughly a sphere, it’s often useful to approximate it by a sphere. So if you’re going to do that, what radius do you use? More generally, what radius do you use when approximating any ellipsoid by a sphere?

This post will discuss the more general problem of finding the radius when approximating any ellipsoid by a sphere. We will give the answer for Earth in particular, and we’ll show how to carry out the calculations. Most of the calculations are easy, but some involve elliptic integrals and we show how to compute these in Python.

## Ellipsoids and spheroids

First of all, what is an ellipsoid? It is a surface whose (x, y, z) coordinates satisfy Earth is an oblate spheroid, which means a = b > c. Specifically, a = b = 6,378,137 meters, and c = 6,356,752 meters.

If you wanted to approximate an ellipsoid by a sphere, you could use

r = (a + b + c)/3.

Why? Because the knee-jerk reaction whenever you need to reduce a set of numbers to one number is to average them.

## Volume of an ellipsoid

We could do a little better, depending on what property of the ellipsoid we’d like to preserve in our approximation. For example, we might want to create a sphere with the same volume as the ellipsoid. In that case we’d use the geometric mean

r = (abc)1/3.

This is because the volume of an ellipsoid is 4πabc/3 and the volume of a sphere is 4πr3/3.

For the particular case of the earth, we’d use

(a2c)1/3 = 6371000.7

## Surface area of an ellipsoid

For some applications we might want a sphere with the same surface area as the ellipsoid rather than the same volume.

The surface area of an ellipsoid is considerably more complicated than the volume. For the special case of an oblate spheroid, like earth, the area is given by where The surface area of a sphere is 4 πr2 and so the following code computes r.

    from math import sqrt, atanh
e = sqrt(1 - (c/a)**2)
r = a*sqrt(1 + (1-e**2)*atanh(e)/e) / sqrt(2)


This gives r = 6371007.1 for the earth, about 6.4 meters more than the number we got matching volume rather than area.

For a general ellipsoid, the surface area is given by where and Here F is the “incomplete elliptic integral of the first kind” and E is the “incomplete elliptic integral of the second kind.” The names are historical artifacts, but the “elliptic” part of name comes from the fact that these functions were discovered in the context of arc lengths with ellipses, so it shouldn’t be too surprising to see them here.

## Computing ellipsoid surface area in Python

In SciPy, F(φ, k) is given by ellipkinc and E(φ, k) is given by ellipeinc. Both function names start with ellip because they are elliptic functions, and end in inc because they are “incomplete.” In the middle, ellipeinc has an “e” because it computes the mathematical function E(φ, k).

But why does ellipkinc have a “k” in the middle? The “complete” elliptic integral of the first kind is K(k) = F(π/2, k). The “k” in the function name is a reminder that we’re computing the incomplete version of the K function.

Here’s the code for computing the surface area of a general ellipsoid:

    from math import sin, cos, acos, sqrt, pi
from scipy.special import ellipkinc, ellipeinc

def area(a, b, c):
phi = acos(c/a)
k = a*sqrt(b**2 - c**2)/(b*sqrt(a**2 - c**2))
E = ellipeinc(phi, k)
F = ellipkinc(phi, k)
elliptic = E*sin(phi)**2 + F*cos(phi)**2
return 2.0*pi*c**2 + 2*pi*a*b*elliptic/sin(phi)


The differences between the various approximation radii are small for Earth. See my next post on elliptical galaxies where the differences are much larger.

Related posts:

# Poor Mercator

This morning’s xkcd cartoon is “What your favorite map projection says about you.” It’s really funny, but poor Mercator comes off as the most boring projection. Mercator is the most familiar projection, but it has some interesting properties. The most important is that lines of constant bearing on the Earth correspond to straight lines on the map, obviously a desirable property for navigation. More details here.

The inverse of the Mercator projection, going from maps onto the globe, is more interesting. Aside from its geographical importance, it gives a way of relating circular and hyperbolic functions without using complex numbers. More details here.

The Mercator projection is also historically interesting. The modern derivation of the Mercator projection uses logarithms and calculus, but Mercator came up with his projection in 1569 before logarithms or calculus had been discovered. (Update: See more historical detail in Thony C’s comment below.)