Spherical coordinate Rosetta Stone

If you’ve only seen one definition of spherical coordinates, you may be shocked to discover that there are multiple conventions. In particular, mathematicians and geoscientists have different conventions. As Volker Michel put it in book on constructive approximation,

Many mathematicians have faced weird jigsaw puzzles with misplaced continents after using a data set from a geoscientist.

Nor are there only two conventions. This post will cover three conventions, which I will call physics, math, and geography. This is not to imply that everyone in a profession uses the conventions named after the profession.

As with Fourier transforms there are many conventions. Presumably the reason for the proliferation of conventions is the same: a useful basic concept will be discovered multiple times in different areas of application. When there are trade-offs between conventions, people working in each area will choose the convention that is simplest for the things they care most about.

Questions and definitions

Given a triple of numbers (c1, c2, c3), that represent spherical coordinates, there are two basic questions:

  1. What do the numbers mean?
  2. What symbol should be used to represent them?

Part of the meaning question is whether the coordinate system is left-handed or right-handed. There is also a minor question of the standard range of equivalent angles.

Polar angle is the angle down from the positive z-axis, the north pole. The north pole has polar angle 0 and the south pole has polar angle π.

Latitude is the angle up from the xy plane or the equator. The latitude of the north pole is π/2 and the latitude of the south pole is −π/2.

Azimuth is the angle from an initial meridian (e.g. the Prime Meridian in geography). In mathematical terms it is the angle between the x-axis and a line connecting the projection of a point to the xy to the origin.

Physics convention

The physics convention is also the ISO 80000-2 convention. The three coordinates are

(radial distance to origin, polar angle, azimuthal angle)

and are conventionally denoted

(r, θ, ϕ).

The point (r, θ, ϕ) corresponds to

(r sin θ cos ϕ, r sin θ sin ϕ, r cos θ)

in Cartesian coordinates.

The coordinate system is right-handed: ϕ increases in the counterclockwise direction when looking down on the xy plane.

The polar angle typically runs from 0 to π and the azimuthal angle typically runs from 0 to 2π.

Math convention

In this convention the three coordinates are

(radial distance to origin, azimuthal angle, polar angle)

and are denoted

(ρ, θ, ϕ).

When compared to the physics notation we have something of a Greek letter paradox. Both conventions agree that the second coordinate should be written θ and the third coordinate should be written ϕ, but the two conventions mean opposite things by the two symbols.

The advantage of using θ to denote the azimuthal angle is that then θ has the same meaning as it has in polar coordinates. Similarly, the advantage of using ρ rather than r is to preserve the meaning of r in polar and cylindrical coordinates: ρ is the distance in 3 dimensions from a point to the origin, and r is the distance in 2 dimensions from the projection of a point to the xy plane and the origin.

Sometimes you’ll see r rather than ρ for the first coordinate.

The point (ρ, θ, ϕ) corresponds to

(ρ cos θ sin ϕ, ρ sin θ sin ϕ, ρ cos ϕ).

The math convention uses the same orientation and typical range coordinate ranges as the physics convention.

Geography and geoscience

Geographic coordinate systems are complicated by the fact that the earth is not perfectly spherical. Geographic latitude and geocentric latitude are not exactly the same as the following highly exaggerated diagram shows.

However, for this post we will assume a perfectly spherical earth in which geographic and geocentric latitudes are the same.

A point on the earth’s surface is described by two numbers, latitude and longitude. This corresponds to spherical coordinates in which the first coordinate is implicitly the earth’s radius. You’ll see latitude and longitude listed in either order.

The polar angle is sometimes called colatitude because latitude is the complementary angle of latitude. Longitude corresponds to azimuth. There are multiple symbols used for latitude, though λ is common for longitude. If you see one coordinate denoted λ, the other is latitude.

The safest, most explicit notation is to use the words longitude and latitude, not depending on symbols or order. Longitude is implicitly longitude east of the Prime Meridian, but to be absolutely clear, put an E at the end. Latitude typically ranges from −π/2 to π/2 (−90° to 90°). Longitude either runs from 0 to 2π or from −π to π (−180° to 180°).

To convert to Cartesian coordinates, set polar angle to π/2 − latitude, azimuth to longitude (east), and use the formulas above for either the physics or mathematics notation. If you have elevation, this is the radial distance above the earth’s surface, so ρ is the earth’s radius plus elevation.

Related posts

How faithful can a map be?

It’s well known that you cannot map a sphere onto the plane without distortion. You can’t map the entire sphere to the plane at all because a sphere and a plane are not topologically equivalent.

But even if you want to map a relatively small portion of globe to paper, say France, with about 0.1% of the globe’s area, you’ll have to live with some distortion. How can we quantify this distortion? How small can it be?

We know the result has to vary with area. We expect it to approach zero as area on the globe goes to zero: the world is not flat, but it is locally flat. And we expect the amount of distortion to go to infinity as we take in more of the globe’s area. We expect to be able to make a map of France without too much distortion, less distortion than making a map of Australia, but more distortion than making a map of Lichtenstein.

This problem was solved a long time ago, and John Milnor wrote an expository article about it in 1969 [1].

Defining distortion

The way to quantify map distortion is to look at the ratio of distances on the globe to distances on paper. We measure the distance between points on the globe by geodesic distance, the length of the shortest path between two points. We’d like the ratio of distances on a globe to distances on a map to be constant, but this isn’t possible. So we look at the minimum and maximum of this ratio. In a good map these two ratios are roughly the same.

Milnor uses

dS(x, y)

to denote the distance between two points x and y on a sphere, and

dE(f(x), f(y))

for the distance between their image under a projection to the Euclidean plane.

The scale of a map with respect to points x and y is the ratio

dE(f(x), f(y)) / dS(x, y).

Let σ1 be the minimum scale as x and y vary and let σ2 be the maximum scale. The distortion of the projection f is defined to be

δ = log(σ21).

When σ1 and σ2 are approximately equal, δ is near 0.

Example: Mercator projection

One feature of the Mercator projection is that there is no distortion along lines of constant latitude. Given two points on the equator (less than 180° apart) their distance on a Mercator projection map is strictly proportional to their distance on the globe. The same is true for two points on the flat part of the US-Canada border along the 49th parallel, or two points along the 38th parallel dividing North Korea and South Korea.

For points along a parallel (i.e. a line of latitude) the scale is constant, such as a thousand miles on the globe corresponding to an inch on the map. We have σ1 = σ2 and so δ = 0.

The distortion in the Mercator projection is all in the vertical direction; distances along a meridian are distorted. Mercator maps latitude φ to something proportional to

φ ↦ log( sec φ + tan φ ).

Near the equator sec φ ≈ 1, tan φ ≈ φ, and the image of φ is approximately φ. But as φ approaches 90° the secant and tangent terms become unbounded and the distortion goes to infinity. You could use the equation for the Mercator projection of latitude to calculate the distortion of a “rectangle on a globe.”

Minimizing distortion

Let Dα be a disk geodesic radius α radians where 0 < α < π. Then Milnor shows that the minimum possible distortion when mapping Dα to the plane is

δ0 = log(α / sin α).

This map is realizable, and is unique up to similarity transformations of the plane. It is known in cartography as the azimuthal equidistant projection.

For small values of α the distortion is approximately α²/6, or 2/3 of the ratio of the area of Dα to the area of the globe.

We said earlier that France takes up about 0.1% of the earth’s surface, and so the minimum distortion for a map of France would be on the order of 0.0007.

[1] John Milnor. A Problem in Cartography. American Mathematical Monthly, December 1969, pp. 1102–1112.

Dutton’s Navigation and Piloting

This morning Eric Berger posted a clip from The Hunt for Red October as a meme, and that made me think about the movie.

I watched Red October this evening, for the first time since around the time it came out in 1990, and was surprised by a detail in one of the scenes. I recognized one of the books: Dutton’s Navigation and Piloting.

Screen shot with Dutton's Navigation and Piloting

I have a copy of that book, the 14th edition. The spine looks exactly the same. The first printing was in 1985, and I have have the second printing from 1989. So it is probably the same edition and maybe even the same printing as in the movie. I bought the book last year because it was recommended for something I was working on. Apparently it’s quite a classic since someone thought that adding a copy in the background would help make a realistic set for a submarine.

My copy has a gold sticker inside, indicating that the book came from Fred L. Woods Nautical Supplies, though I bought my copy used from Alibris.

Here’s a clip from the movie featuring Dutton’s.

Dutton’s has a long history. From the preface:

Since the first edition of Navigation and Nautical Astronomy (as it was then titled) was written by Commander Benjamin Dutton, U. S. Navy, and published in 1926, this book has been updated and revised. The title was changed after his death to more accurately reflect its focus …

The 14th edition contains a mixture of classical and electronic navigation, navigating by stars and by satellites. It does not mention GPS; that is included in the latest edition, the 15th edition published in 2003.

Related posts

Area of a “rectangle” on a globe

What do we mean by rectangle?

This post will derive the area of a spherical region bounded by lines of latitude and longitude. Such a region corresponds to an actual rectangle on a Mercator projection map, with sides aligned with the coordinate axes, and is approximately a rectangle on a sphere if the rectangle is not too big [1].

What do we know up front?

Before we get into detailed equations, we know that the area will be proportional to the difference in longitude. If we’re looking that the area between two parallels, such as the equator and the Arctic Circle, the area between 10° and 20° longitude is the same as the area between 80° and 90° longitude, and twice the area between 72° and 77° longitude.

The difficulty is latitude. Say we look at squares on a map that are 1° of longitude wide and 1° of latitude tall. Those squares are on the map will correspond to more area on the globe for latitudes near the equator, and less area at high latitudes.

So the area bounded by longitudes θ1 and θ2 and latitudes φ1 and φ2 will depend on φ1 and φ2 individually, but only on the difference θ1 − θ2.

Spherical caps

The region on a sphere above a fixed line of latitude is called a spherical cap. The northern hemisphere, the region above the equator, would be a very large spherical cap. The region inside the Arctic Circle would be a smaller spherical cap.

Let R be the radius of the earth. Then the surface area above a latitude φ is

A = 2πR²(1 − sin φ).

You could derive this using calculus by thinking of the spherical cap as a surface of revolution.

Spherical bands

Given two latitudes φ1 and φ2 with latitudes φ1 > φ2, the area of a band between latitude φ1 and latitude φ2 is the area of the spherical cap above latitudes φ2 minus the area of the spherical cap above latitudes φ1. This gives

A = 2πR²(sin φ1 − sin φ2).

Area of latitude/longitude grid

Now we can find the area of the region bounded by longitudes θ1 and θ2 and latitudes φ1 and φ2. The total area between latitudes φ1 and φ2 is given by the equation above. The subset of this area between longitudes θ1 and θ2 is proportional to θ1 − θ2. If longitude is measured in radians then

A = R² (sin φ1 − sin φ2) (θ1 − θ2).

If longitude is measured in degrees, we have

A = π R² (sin φ1 − sin φ2) (θ1 − θ2)/180.

Related posts

[1] As Nathan helpfully pointed out in the comments, it might be clearer to say “not too big and not too close to one of the poles.” Another way to put it might be to say the idea of what is “big” depends on how close to a pole you are.

Law of cosines on a sphere

The previous post looked at the analog of the Pythagorean theorem on a sphere. This post looks at the law of cosines on a sphere.

Yesterday we looked at triangles on a sphere with sides a and b meeting at a right angle and hypotenuse c. Denote the angle opposite a side with the capital version of the same letter, so C is the angle opposite c. We assumed C is a right angle, but now we will remove that assumption.

The spherical analog of the law of cosines says

cos(c) = cos(a) cos(b) + sin(a) sin(b) cos(C).

Note that we have two kinds of angles: the arcs that make up the sides, and the angle formed by the intersection of arcs [1]. So the cos(C) term at the end is different animal from the other terms in the equation.

If C is a right angle, cos(C) = 0 and the second term drops out, leaving us with the spherical counterpart of the Pythagorean theorem. But we do not require that C be a right angle.

Application to air distance

Suppose we want to find how far a plane would travel between Los Angeles (LAX) and Houston (IAH), assuming it takes a great circle path. The lat/long coordinates of the two airports are (33.94°, -118.41°) and (29.98°, -95.34°).

Los Angeles and Houston are approximately at the same latitude, but even if they were at exactly the same latitude we couldn’t just find the flight distance by finding the length of the arc of constant latitude between them because that arc would not be part of a great circle.

To find the distance between LAX and IAH we form a triangle with vertices at LAX, IAH, and the north pole. Call the arc from LAX to the north pole a and the arc from IAH to the north pole b. Since latitude is the angle up from the equator, the angle down from the pole is 90° minus latitude. So

a = 90° – 33.94° = 56.06°

and

b = 90° – 29.98° = 60.02°

The arcs a and b meet at the angle C equal to the differences in the two longitudes. That is,

C = 118.41° – 95.34° = 23.07°.

The law of cosines

cos(c) = cos(a) cos(b) + sin(a) sin(b) cos(C)

reduces to

cos(c) = cos(56.06°) cos(60.02°) + sin(56.06°) sin(60.02°) cos(23.07°)

in our problem.

This tells us

cos(c) = 0.9402

from which we find c = 0.3477 radians or 19.92°. Assume the earth is a perfect sphere of radius r = 3959 miles. The arc length we’re after is r times c in radians, or 1376 miles.

Related posts

[1] In the language of differential geometry, the arcs are geodesics on the sphere and the angles of intersection are measured in the tangent space at the intersection point.

Pythagoras on a sphere

Suppose you drive a distance a, turn right, then drive a distance b. How far are you from where you started?

If a and b are relatively small, then the answer is given by the Pythagorean theorem. Small here means small relative to the radius of the Earth. Any distance you drive is small enough: it matters more that roads are not perfectly straight than that the Earth is not perfectly flat.

Now suppose you fly a distance a, turn right, then fly a distance b. Now how far are you from where you started? Now the curvature of the earth might matter a good deal, and the Pythagorean theorem could give you a poor answer.

Spherical version of the Pythagorean theorem

If a and b are moderately large, we need to look at spherical triangles, not plane triangles. A “straight” line on a sphere is a piece of a great circle, a circle whose center is the center of the sphere. A spherical triangle has three straight sides, straight in the sense of a straight line on a sphere. We need to consider a spherical triangle with sides a and b that meet at right angles, and we want to find the length of the hypotenuse c of the spherical triangle.

The spherical analog of the Pythagorean theorem says

cos(c) = cos(a) cos(b).

This looks weird the first time you see it, because you’re taking the cosine of a side rather than an angle. But on a sphere, the sides are angles. For example, if a side of the triangle is on a meridian, we could measure its length in degrees of latitude. But the theorem above does not require sides to be lined up with any coordinate system. It only requires that the sides a and b meet at a right angle.

Small spherical triangles

Let’s suppose for a moment that a and b are small when measured in radians, small enough that a², b², and ab are negligible compared to a and b. Then we can use the approximation

cos(x) ≈ 1 – x²/2

in the theorem above to show that

c² ≈ a² + b².

This is reassuring since we know the Pythagorean theorem holds for human-scale distances on our planet.

Example with a big triangle

Quito to Nairobi to Jerusalem and back to Quito

Now let’s look at a big triangle. The vertices have longitude and latitude (-80°, 0°), (35°, 0°), and (35°, 30°). These correspond approximately to Quito, Nairobi, and Jerusalem. So a = 115° and b = 30°. We can solve for c and find it equals 111°.

A minute of latitude is a nautical mile; the nautical mile is defined the way to make these kinds of calculations easier. A nautical mile is about 15% longer than a mile.

A degree of a great circle is 60 nautical miles. So our hypotenuse of 111 degrees is 6,660 nautical miles. If we had used the Pythagorean theorem, we would have overestimated this distance as 7,130 nautical miles.

Notice that the triangle with vertices at Quito, Nairobi, and Jerusalem doesn’t look like a triangle when projected flat. Two legs look straight because they follow lines of constant latitude or longitude. But the flight from Jerusalem to Quito is just as straight as the other flights if you look at it on a globe.

Not so great circles

In the example above, I chose Quito and Nairobi because they were near the equator. The equator is a great circle.

Suppose I’d chosen Houston and Wuhan instead. Both cities are at about 30° latitude. I could not find the distance between the two cities by simply subtracting their longitudes because the 30th parallel is not a great circle. Spherical geometry is generally more involved than in the example above.

A great circle path from Houston to Wuhan would not follow the 30th parallel but would actually cross Alaska. The great circle distance is 6,800 nautical miles. Flying along a path of 30° latitude would add an extra 1,000 nautical miles to the flight.

Related posts

A tale of three cities

Big blue marble

Pick three cities and form a spherical triangle by connecting each pair of cities with the shortest arc between them. How might you find the area of this triangle?

For this post, I’ll assume the earth is perfectly spherical. (Taking into account that the earth is slightly oblate makes the problem much more complicated. Maybe a post for another time.)

Spherical excess

If you draw a triangle on the plane, the interior angles add up to 180°. But if you draw a triangle on a sphere, the interior angles add up to more than 180°. The amount by which the sum of the interior angles exceeds 180° is called the spherical excess. It turns out that the area of a spherical triangle is proportional to its spherical excess.

For example, if a spherical triangle is small relative to the radius of the sphere, it’s nearly a Euclidean triangle, and the spherical excess is negligible. But consider a much bigger spherical triangle, one with a vertex at the north pole and vertices at two points 90° apart on the equator. This triangle has three right angles, so the sum of the interior angles is 270° and the spherical excess is 90°.

Relating spherical excess and area

If we measure spherical excess E in radians, then the area of a spherical triangle is

A = ER²

where R is the radius of the sphere. This remarkable theorem was discovered by Thomas Harriot in 1603.

Let’s go back to the example above of a triangle running from the north pole to two points 90° apart on the equator. This triangle takes up 1/8 of the earth’s surface, and area of the earth is 4πR², and so the triangle has area πR²/2. In this case the spherical excess E was π/2, and so we could have come to the same result by multiplying E by R².

Computing area

We can find the area of a spherical triangle by measuring the angle of arcs between pairs of points. We’ll compute area assuming the sphere has radius 1. (Otherwise temporarily assume radius 1 and then multiply by R² when we’re done.)

The previous post explains how to find a parameterization for the great circle connecting points on a sphere. We can take derivatives to find tangent vectors where the great circles intersect, and compute the angle between these tangents to find the interior angles of our triangles.

The previous post showed that if the central angle between two vectors A and B is θ then

cos(t) A + sin(t) (cot θ A – csc θ B)

parameterizes a great circle including A and B. This curve passes through A at time t = 0 with velocity

cot θ A – csc θ B

and so this gives us a tangent vector at A in the direction of B. Repeat the exercise for a great circle between A and C. Then the cosine of the angle of intersection is the two normalized tangent vectors. We can thus obtain all the interior angles, and from there we can compute the spherical excess and thus the area of the spherical triangle.

Related posts

Maidenhead geocode system

The Maidenhead Locator System encodes a pair of longitude and latitude coordinates in a slightly complicated but ingenious way. Amateur radio operators using this geocoding system to describe locations.

The Wikipedia article on the subject describes the what of the system, but I’d like to say more about the why of the system. I’ll also go through an example in great detail. Continue reading

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)

Research Triangle

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.

More spherical trig posts

Ellipsoid distance on Earth

globe

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