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 data-sets 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:

distance-amsterdam-london.png

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:

map-flat.png

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 anti-meridian at 180°W/E). 5.03° in East-West difference is then about 40 000 ✕ (5.03/360) = 559 km and the North-South 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.

sphere2.png

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 great-circle 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:

great-circle.jpg

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 geo-located 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.

Shortlink

This article has a short URL available: http://drck.me/spat-dist-8kf

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:

// Helsinki-Malmi airport (EFHF)
 $efhf = new midgardmvc_helper_location_spot(60.254558, 25.042828);
 // Helsinki-Vantaa 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).

Cool stuff, thank you!

Looking forward to the MySQL specific things as the web seems to be lacking good articles on it

Great write-up!

On a totally unrelated note, the PHP-code snippets you encapsuled with the &lt;?php and ?&gt; don't show up in the RSS-feed. Normal &lt;code&lt; does, though.

Also, displaying the footnotes even below the comment-section is quite unusual und unexpected.

Nevertheless, great write-up!

This site has a great summary of methods to calculate geographical distances:

http://www.movable-type.co.uk/scripts/latlong.html

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.

Add Comment

Name:
Email:

Will not be posted. Please leave empty instead of filling in garbage though!
Comment:

Please follow the reStructured Text format. Do not use the comment form to report issues in software, use the relevant issue tracker. I will not answer them here.


All comments are moderated
  • 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.

Life Line