Spatial Indexes: Calculating Distance
During my "Geolocation and Mapping with PHP" talk that I've given a few times I briefly touch on the subject of indexes on datasets of spatial data. This isn't as simple as just solving Pythagoras theorem and this article is meant to clarify this.
The flat Earth model
Pythagoras theorem can be used to calculate the distance between two points quite easily; you take the square root of the square of the absolute vertical distance plus the square of the absolute horizontal distance; in short:
d = β(x1  x2Β² + y1  y2Β²)
or in PHP:
$d = sqrt(pow(abs($x1  $x2), 2) + pow(abs($y1  $y2), 2));
If you take for example London's coordinates (51.50Β°N, 0.13Β°W) and Amsterdam's coordinates (52.37Β°N, 4.90Β°E) we can calculate the distance with:
d = β(0.13  4.90Β² + 51.50  52.37Β²) d = β(5.03Β² + 0.87Β²) d = β(26.0578) d β 5.10
And on a map:
But what does a difference of "5.10Β°" actually mean? How far is this in useful units, such as meters?
If we show the whole map of which the above is an extract, we come to:
The distance around the equator, and through the poles is roughly the same, 40.000km (please be aware that the blue line only shows half of it, the other half is going through the antimeridian at 180Β°W/E). 5.03Β° in EastWest difference is then about 40 000 β (5.03/360) = 559 km
and the NorthSouth difference about 20 000 β (0.87/180) = 97 km
. Using those numbers within the Pythagoras theorem we end up with a distance of β(559Β² + 97Β²) = 567 km
. Although the calculation is correct, the answer is still wrong. The real distance is closer to 360 km.
The spherical Earth model
If we look at the Earth in its original (mostly) spherical 1 shape, then it's clear that 10Β° longitude (East/West) at 60Β°N is going to be less of a distance than 10Β°E/W at the equator. It's actually fairly easy to calculate how much 1Β° longitude is at 60Β°N by using cos(60) * 1/360 * 6371 * 2Ο
2. Or in PHP:
<?php $oneDeg = cos(deg2rad(60)) * // adjustment for latitude and radians/degrees 1/360 * // 1 out of 360Β° 6371 * 2 * M_PI, // circumference of the Earth at the equator "\n"; ?>
This returns 55.597 km per Β°
for 60Β°N. The same distance in degrees on the Equator gives 6371 * 2 * M_PI / 360 = 111.195
.
The following diagram shows ones more that latitudinal degrees always correspond with the same distance in kilometer, whereas longitudinal degrees differ.
In the diagram the line A is a line from 0Β°N, 90Β°W to 10Β°N, 90Β°W. It has the same length as line B, from 30Β°N, 90Β°W to 40Β°N, 90Β°W: a 36th of the circumference of the Earth through the poles. Line C, from 0Β°N, 30Β°W to 0Β°N, 20Β°W has the same length. Line D however, from 50Β°N, 30Β°W to 50Β°N, 20Β°W is shorter by a factor of cos(50Β°) β
0.64
.
If we look again at the distance between London and Amsterdam, ignore the differences in latitude and instead pick the average, we see:
<?php $d = abs(0.13  4.90); // $d = 5.03 degrees $e = 5.03/360 * cos(deg2rad(51.935)) * 6371 * 2 * M_PI; $e = 5.03 * 0.617 * 6371 * 2 * M_PI; // $e β 345 km ?>
Which is a bit shorter than the expected 360km, but that's because we conveniently forgot about the difference in latitude.
Sadly, we can't use Pythagoras's theorem to calculate the real distance with the latitude difference taken account as well. This is because the theorem is meant for Euclidean geometry, and a sphere does not follow the rules of this geometry. Instead we need to use a formula that is called the greatcircle distance formula. The main concept behind it is that a circle is drawn across the whole sphere that connects both the start (point P) as well as the end (point V). Then with that circle the distance can be calculated. The following diagram shows this:
I will spare you how the function is derived, but the distance calculation ends up as being:
<?php function distance($latA, $lonA, $latB, $lonB) { // convert from degrees to radians $latA = deg2rad($latA); $lonA = deg2rad($lonA); $latB = deg2rad($latB); $lonB = deg2rad($lonB); // calculate absolute difference for latitude and longitude $dLat = ($latA  $latB); $dLon = ($lonA  $lonB); // do trigonometry magic $d = sin($dLat/2) * sin($dLat/2) + cos($latA) * cos($latB) * sin($dLon/2) *sin($dLon/2); $d = 2 * asin(sqrt($d)); return $d * 6371; } ?>
If we punch in our original numbers form London (51.50Β°N, 0.13Β°W) and Amsterdam (52.37Β°N, 4.90Β°E), we calculate the following:
<?php $d = distance(51.50, 0.13, 52.37, 4.90); echo $d, " km\n"; ?>
Which gives us the expected result of 358.07 km
.
Conclusion
I hope that the above clarified the difference between 2D spatial indexing with the flat Earth model and spatial indexing of geolocated data (the spherical Earth model). In future articles I will go into specific implementations of spatial indexing by traditional databases such as SQLite, MySQL and PostGreSQL; NoSQL databases such as MongoDB and CouchDB; and Solr.
Comments
In the first section, you have d = β(x1  x2Β² + y1  y2Β²)
but why the need for the absolute values? Squaring the difference will result in a positive number no matter what, so shave a smidge off. Unless there's some other fundamental reason for it.
@Glen: You're right; there is absolutely no reason for that :)
Just what I looked for right now. Thanks a lot. Saved me some time today.
Btw, here is the same code in python:
def getDistance(pointA, pointB): # convert from degrees to radians latA = math.radians(pointA[0]) lonA = math.radians(pointA[1]) latB = math.radians(pointB[0]) lonB = math.radians(pointB[1]) #calculate absolute difference for latitude and longitude dLat = (latA  latB) dLon = (lonA  lonB) # do trigonometry magic d = math.sin(dLat/2) * math.sin(dLat/2) + math.cos(latA) \ * math.cos(latB) * math.sin(dLon/2) * math.sin(dLon/2) d = 2 * math.asin(math.sqrt(d)) return (d * 6371)
...or you can just:
// HelsinkiMalmi airport (EFHF) $efhf = new midgardmvc_helper_location_spot(60.254558, 25.042828); // HelsinkiVantaa airport (EFHK) $efhk = new midgardmvc_helper_location_spot(60.317222, 24.963333); // There are 8.2 kilometers between the airports $distance = $efhf>distance_to($efhk); echo $distance; // 8.2
(I know, I should really move this library to Zeta Components).
Another Python implementation:
https://github.com/cannonerd/adventure_tablet/blob/master/src/point.py
Cool stuff, thank you!
Looking forward to the MySQL specific things as the web seems to be lacking good articles on it
Great writeup!
On a totally unrelated note, the PHPcode snippets you encapsuled with the <?php and ?> don't show up in the RSSfeed. Normal <code< does, though.
Also, displaying the footnotes even below the commentsection is quite unusual und unexpected.
Nevertheless, great writeup!
This site has a great summary of methods to calculate geographical distances:
Interesting post. We spoke about this on the drupal dev days in belgium: http://bxl2011.drupaldays.org/node/304
I am very interested in your follow up for solr! In drupal there is geo search for mysql. But this distance as a solr result is not yet possible (out of the box).
Keep us informed!
Great !
Excellent, clear, concise, and very useful.

1
The Earth is not really a sphere, but an approximation of it. However, doing the same calculations for an ellipsoid can (as far as I know) only be done by approximation. The difference would hardly matter for finding the closest pub.

2
In this article, I've used an average radius of the Earth of 6371 km.
Shortlink
This article has a short URL available: https://drck.me/spatdist8kf