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
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.
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%.
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
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:
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.
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.
What point on Earth is farthest from its center? Mt. Everest comes to mind, but its summit is the point highest above sea level, not the point farthest from the center. These are not the same because the Earth is not perfectly spherical.
Our planet bulges slightly at the equator due to its rotation. The equatorial diameter is about 0.3% larger than the polar diameter. Sea level at the equator is about 20 kilometers farther from the center of the Earth than sea level at the poles.
Mt. Everest is about nine kilometers above sea level and located about 28 degrees north of the equator. Chimborazo, the highest point in Ecuador, is seven kilometers above sea level and 1.5 degrees south of the equator.
So how far are Mt. Everest and Chimborazo from the center of the Earth? To answer that, we first need to how far sea level at latitude θ is from the center of the Earth.
Imagine slicing the Earth with a plane containing its polar diameter. To a good approximation (within 100 meters) the resulting shape would be an ellipse. The equation of this ellipse is
(x / a)2 + (y / b)2 = 1
where a = 6378.1 km is the equatorial radius and b = 6356.8 km is the polar radius. A line from the center of the ellipse to a point at latitude θ has equation y = tan(θ) x. Solving the pair of equations for x shows that the distance from the center to the point at latitude θ is
d = sqrt( a2b2 sec2 θ / (a2 tan2 θ + b2 ) )
For Mt. Everest, θ = 27.99 degrees and so d = 6373.4. For Chimborazo, θ = -1.486 degrees and so d = 6378.1. So sea level is 4.7 km higher at Chimborazo. Mt. Everest is 2.6 km taller, but the summit of Chimborazo is about 2.1 km farther away from the center of the Earth.
Update: See my next post for a slight correction. A more accurate calculation would compute sea level is about 4.65 km higher at Chimborazo than Mt. Everest.
Universal time (UTC) is the same as Greenwich Mean Time (GMT), give or take a second. It’s essentially the time in Greenwich, England except it ignores Daylight Savings Time.
The abbreviation UTC is an odd compromise. The French wanted to use the abbreviation TUC (temps universel coordonné) and the English wanted to use CUT (coordinated universal time). The compromise was UTC, which doesn’t actually abbreviate anything.
Sometimes a ‘Z’ is appended to a time to indicate it is expressed in UTC. The NATO phonetic alphabet code for ‘Z’ is ZULU, and so UTC is sometimes called “Zulu Time.”
It’s useful to know how your time zone relates to UTC. (You can look it up here.) For example, I live in the US Central time zone. Central Standard Time (CST) is UTC-6, i.e. we’re 6 hours behind Greenwich. Knowing your time relative to UTC makes international communication easier. It also helps you read computer time stamps since these almost always use UTC.
One of the advantages of UTC is that it avoids Daylight Savings Time. DST is surprisingly complicated when you look at it in detail. Some countries observe DST and some do not. Countries that do observe DST may begin and end DST on different dates, and those dates can change from year to year. And inside countries that observe DST some regions are exceptions. For example, the United States generally observes DST, but Arizona does not. Actually, it’s even more complicated: The Navajo Nation inside Arizona does observe DST.
A natural approach to mapping the Earth is to imagine a cylinder wrapped around the equator. Points on the Earth are mapped to points on the cylinder. Then split the cylinder so that it lies flat. There are several ways to do this, all known as cylindrical projections.
One way to make a cylindrical projection is to draw a lines from the center of the Earth through each point on the surface. Each point on the surface is then mapped to the place where the line intersects the cylinder. Another approach would be to make horizontal projections, mapping each point on Earth to the closest point on the cylinder. The Mercator projection is yet another approach.
With any cylindrical projection parallels, lines of constant latitude, become horizontal lines on the map. Meridians, lines of constant longitude, become vertical lines on the map. Cylindrical projections differ in how the horizontal lines are spaced. Different projections are useful for different purposes. Mercator projection is designed so that lines of constant bearing on the Earth correspond to straight lines on the map. For example, the course of a ship sailing northeast is a straight line on the map. (Any cylindrical projection will represent a due north or due east course as a straight line, but only the Mercator projection represents intermediate bearings as straight lines.) Clearly a navigator would find Mercator’s map indispensable.
Latitude lines become increasingly far apart as you move toward the north or south pole on maps drawn with the Mercator projection. This is because the distances between latitude lines has to change to keep bearing lines straight. Mathematical details follow.
Think of two meridians running around the earth. The distance between these two meridians along a due east line depends on the latitude. The distance is greatest at the equator and becomes zero at the poles. In fact, the distance is proportional to cos(φ) where φ is the latitude. Since meridians correspond to straight lines on a map, east-west distances on the Earth are stretched by a factor of 1/cos(φ) = sec(φ) on the map.
Suppose you have a map that shows the real time position of a ship sailing east at some constant rate. The corresponding rate of change on the map is proportional to sec(φ). In order for lines of constant bearing to be straight on the map, the rate of change should also be proportional to sec(φ) as the ship sails north. That says the spacing between latitude lines has to change according to h(φ) where h‘(φ) = sec(φ). This means that h(φ) is the integral of sec(φ) which equals log |sec(φ) + tan(φ)|. The function h(φ) becomes unbounded as φ approaches ± 90°. This explains why the north and south poles are infinitely far away on a Mercator projection map and why the area of northern countries is exaggerated.
(Update: The inverse of the function h(φ) has some surprising properties. See Inverse Mercator projection.)
The modern explanation of Mercator’s projection uses logarithms and calculus, but Mercator came up with his projection in 1569 before logarithms or calculus had been discovered.
For more details of the Mercator projection, see Portraits of the Earth.
Someone sent me email regarding my online calculator for computing the distance between to locations given their longitude and latitude values. He wants to do sort of the opposite. Starting with the longitude and latitude of one location, he wants to find the longitude and latitude of locations moving north/south or east/west of that location. I like to answer reader questions when I can, so here goes. I’ll give a theoretical derivation followed by some Python code.
Longitude and latitude and are usually measured in degrees, but theoretical calculations are cleaner in radians. Someone using the Python code below could think in terms of degrees; radians will only be used inside function implementations. We’ll use the fact that on a circle of radius r, an arc of angle θ radians has length rθ. We’ll assume the earth is a perfect sphere. See this post for a discussion of how close the earth is to being a sphere.
I’ll start with moving north/south since that’s simpler. Let R be the radius of the earth. An arc of angle φ radians on the surface of the earth has length M = Rφ, so an arc M miles long corresponds to an angle of φ = M/R radians. Moving due north or due south does not change longitude.
Moving east/west is a little more complicated. At the equator, the calculation is just like the calculation above, except that longitude changes rather than latitude. But the distance corresponding to one degree of longitude changes with latitude. For example, one degree of longitude along the Arctic Circle doesn’t take you nearly as far as it does at the equator.
Suppose you’re at latitude φ degrees north of the equator. The circumference of a circle at constant latitude φ, a circle parallel to the equator, is cos φ times smaller than the circumference of the equator. So at latitude φ an angle of θ radians describes an arc of length M = R θ cos φ. A distance M miles east or west corresponds to a change in longitude of θ = M/(R cos φ). Moving due east or due west does not change latitude.
The derivation above works with angles in radians. Python’s cosine function also works in radians. But longitude and latitude are usually expressed in degrees, so function inputs and outputs are in degrees.
import math
# Distances are measured in miles.
# Longitudes and latitudes are measured in degrees.
# Earth is assumed to be perfectly spherical.
earth_radius = 3960.0
degrees_to_radians = math.pi/180.0
radians_to_degrees = 180.0/math.pi
def change_in_latitude(miles):
"Given a distance north, return the change in latitude."
return (miles/earth_radius)*radians_to_degrees
def change_in_longitude(latitude, miles):
"Given a latitude and a distance west, return the change in longitude."
# Find the radius of a circle around the earth at given latitude.
r = earth_radius*math.cos(latitude*degrees_to_radians)
return (miles/r)*radians_to_degrees
I posted a page this evening that lets you calculate the distance between two locations using their latitudes and longitudes. I’ve had to do this calculation once in a while and thought I’d make it available online for anyone else who needed to do the same. There is one page providing an online calculator and another page giving the formula used to calculate the distances and its derivation.
This project was prompted by a friend asking me how far my home is from Galveston where the hurricane is supposed to make landfall this weekend. Since we live northwest of Houston, we’re pretty far inland.