PHP – Calculating Distance Between Two Locations Given Their GPS Coordinates

Amongst other things, I have the habit of geotagging my photographs and for this purpose, I use my Garmin GPSMap 60CSx, an API I wrote to store tracks in MySQL & Geonames data in a MySQL [please refer to this if you wish to do so as well]. For completeness sake, the MySQL table structure that holds Geonames data is given below:

DROP TABLE IF EXISTS `MyDatabase`.`allCountries` ;
CREATE TABLE IF NOT EXISTS `MyDatabase`.`allCountries` (
  `geo_id`               INT(11) UNSIGNED NOT NULL PRIMARY KEY,
  `geo_name`             VARCHAR(200) NOT NULL DEFAULT '',
  `geo_ansiname`         VARCHAR(200) NOT NULL DEFAULT '',
  `geo_alternate_names`  VARCHAR(2000) NOT NULL DEFAULT '',
  `geo_latitude`         DOUBLE PRECISION(11,7) NOT NULL DEFAULT '0',
  `geo_longitude`        DOUBLE PRECISION(11,7) NOT NULL DEFAULT '0',
  `geo_feature_class`    CHAR(1) ,
  `geo_feature_code`     VARCHAR(10) ,
  `geo_country_code`     CHAR(2),
  `geo_country_code2`    VARCHAR(60),
  `geo_admin1_code`      VARCHAR(20) DEFAULT '',
  `geo_admin2_code`      VARCHAR(80) DEFAULT '',
  `geo_admin3_code`      VARCHAR(20) DEFAULT '',
  `geo_admin4_code`      VARCHAR(20) DEFAULT '',
  `geo_population`       BIGINT(11) DEFAULT '0',
  `geo_elevation`        INT(11) DEFAULT '0',
  `geo_gtopo30`          INT(11) DEFAULT '0',
  `geo_timezone`         VARCHAR(40),
  `geo_mod_date`         DATE DEFAULT '0000-00-00'

Underlying Mathematics

Something which most of us learn in high school mathematics classes is the well known fact that the shortest distance between any two given points is a straight line. If \vec{P}\:\left(x_1, y_1, z_1\right) and \vec{Q}\:\left(x_2, y_2, z_2\right) are the Cartesian coordinates of two points, the distance between them, d, would be:

d \;=\; \sqrt{ \left(x_2 \:-\: x_1 \right)^2 \:+\: \left(y_2 \:-\: y_1 \right)^2 \:+\: \left(z_2 \:-\: z_1 \right)^2}

Spherical Polar CoordinatesWhile dealing with Cartesian coordinates looks plain and simple, it’s not all that helpful for our intended purpose – for what it’s worth, we live on a curved surface and for this, Spherical Polar coordinates come in really handy. And the shortest distance between any two points on the surface of a sphere is no longer a straight line but a great circle [read: geodesic]. If r is the radius of our beautiful Earth, \theta the elevation angle (almost similar to latitude – almost, yeah; explanation comes soon) and \phi the azimuthal angle (similar to longitude), then we can write the points \vec{P}\:\left(x_1, y_1, z_1\right) and \vec{Q}\:\left(x_2, y_2, z_2\right) as \vec{P}\:\left(r, \theta_1, \phi_1\right) and \vec{Q}\:\left(r, \theta_2, \phi_2\right). The transformation from Cartesian to Spherical Polar coordinates is not that difficult and may be written as

x_1 \;=\; r\:\sin\theta_1\:\cos\phi_1 \hspace{0.50in} y_1 \;=\; r\:\sin\theta_1\:\sin\phi_1 \hspace{0.50in} z_1 \;=\; r\:\sin\theta_1

x_2 \;=\; r\:\sin\theta_2\:\cos\phi_2 \hspace{0.50in} y_2 \;=\; r\:\sin\theta_2\:\sin\phi_2 \hspace{0.50in} z_2 \;=\; r\:\sin\theta_2

If the straight lines drawn through these points, \vec{P}\:\left(r, \theta_1, \phi_1\right) and \vec{Q}\:\left(r, \theta_2, \phi_2\right), intersect each other at the center of the Earth and subtend an angle \alpha, then by extending the Law of Cosines [from plane trigonometry] to the spherical case, it can be shown that

\cos\alpha \;=\; \left(1/r^2\right) \left( \vec{P} \:\cdot\: \vec{Q} \right)

After some simplification and using basic trigonometric identities, one can arrive at

\cos\alpha \;=\; \cos\theta_1 \cdot \cos\theta_2 \:+\: \sin\theta_1 \cdot \sin\theta_2 \cdot \cos\left(\phi_1 - \phi_2\right)

It is generally assumed that elevation angle, \theta, can take values in the range \left[ 0, \pi\right] and azimuthal angle, \phi, can take values in the range \left[ 0, 2\pi\right]. However, the accepted geographic convention says that latitude (elevation angle) can range between \left[ -\pi/2, \pi/2\right] and that longitude (azimuth angle) can range between \left[ -\pi, \pi\right]. The expression we have above for \cos\alpha is still valid, except that we need to replace \theta with \left(\pi/2 - \theta\right) and \phi by \left(\pi - \phi\right). Once those replacements are made, the result becomes

\cos\alpha \;=\; \sin\theta_1 \cdot \sin\theta_2 \:+\: \cos\theta_1 \cdot \cos\theta_2 \cdot \cos\left(\phi_1 - \phi_2\right)

\alpha \;=\; \arccos \left[ \sin\theta_1 \cdot \sin\theta_2 \:+\: \cos\theta_1 \cdot \cos\theta_2 \cdot \cos\left(\phi_1 - \phi_2\right) \right]

where \theta indeed represents the latitude and \phi the longitude. It is a well known fact that there are 60 nautical miles per degree of latitude [number of nautical miles per degree of longitude depends on the value of latitude at that place] and that there are 1.1515 statute (or regular) miles in every nautical mile. If \alpha is expressed in degrees (not radians), then the distance between points \vec{P}\:\left(r, \theta_1, \phi_1\right) and \vec{Q}\:\left(r, \theta_2, \phi_2\right) will be

d \;=\; \arccos \left[ \sin\theta_1 \cdot \sin\theta_2 \:+\: \cos\theta_1 \cdot \cos\theta_2 \cdot \cos\left(\phi_1 - \phi_2\right) \right] \times 60 \mbox{\hspace{0.10in}nautical miles}

d \;=\; \arccos \left[ \sin\theta_1 \cdot \sin\theta_2 \:+\: \cos\theta_1 \cdot \cos\theta_2 \cdot \cos\left(\phi_1 - \phi_2\right) \right] \times 60 \times 1.1515 \mbox{\hspace{0.10in}miles}

This looks nice, elegant and easily computable but has an inherent problem: since \arccos() is not so well behaved, the resulting distance could be erroneous if the longitudinal separation between \vec{P}\:\left(r, \theta_1, \phi_1\right) and \vec{Q}\:\left(r, \theta_2, \phi_2\right) is not large enough. As such, it can become a significant problem when my API is trying to look up the database in search of nearest places – for geotagging purposes.

A slightly more complicated looking approach/formula that solves this issue is the haversine formula. One can start from Spherical Law of Cosines and use the standard trigonometric identities,

\cos\theta \:=\: 1 - \sin^2\left(\theta/2\right) \hspace{0.25in} \mbox{and} \hspace{0.25in} \cos\left(\alpha + \beta\right) \:=\: \cos\alpha\:\cos\beta \:-\: \sin\alpha\:\sin\beta

and arrive at the Haversine Formula (R being the radius of our very own earth):

\mbox{haversin}\left(d/R\right) \;=\; \mbox{haversin}\:\Delta\theta \:+\: \cos\theta_1 \cdot \cos\theta_2 \cdot \mbox{haversin}\:\Delta\phi

\sin^2\left(d/2R\right) \;=\; \sin^2\left(\Delta\theta/2\right) \:+\: \cos\theta_1 \cdot \cos\theta_2 \cdot \sin^2\left(\Delta\phi/2\right)

\sin\left(d/2R\right) \;=\; \sqrt{\sin^2\left(\Delta\theta/2\right) \:+\: \cos\theta_1 \cdot \cos\theta_2 \cdot \sin^2\left(\Delta\phi/2\right)}

d \;=\; 2R\cdot \arcsin\left(\sqrt{\sin^2\left(\Delta\theta/2\right) \:+\: \cos\theta_1 \cdot \cos\theta_2 \cdot \sin^2\left(\Delta\phi/2\right)}\right)

d \;=\; 2R\cdot \arcsin\left[\min\left(1, \sqrt{\sin^2\left(\Delta\theta/2\right) \:+\: \cos\theta_1 \cdot \cos\theta_2 \cdot \sin^2\left(\Delta\phi/2\right)}\right)\right]

The \min() in last step is a result of the fact that \sin() maxes out at 1 and there is a non-zero probability that the square root thing can exceed 1. While not necessarily perfect – especially when \vec{P}\:\left(r, \theta_1, \phi_1\right) and \vec{Q}\:\left(r, \theta_2, \phi_2\right) are antipodal – this does provide a pretty accurate estimate for distance.

PHP – Spherical Law of Cosines

<!--?php $earth_radius = 3960.00; # in miles $lat_1 = "47.117828"; $lon_1 = "-88.545625"; $lat_2 = "47.122223"; $lon_2 = "-88.568781"; $delta_lat = $lat_2 - $lat_1 ; $delta_lon = $lon_2 - $lon_1 ; # Spherical Law of Cosines function distance_slc($lat1, $lon1, $lat2, $lon2) { global $earth_radius; global $delta_lat; global $delta_lon; $distance = sin(deg2rad($lat1)) * sin(deg2rad($lat2)) + cos(deg2rad($lat1)) * cos(deg2rad($lat2)) * cos(deg2rad($delta_lon)) ; $distance = acos($distance); $distance = rad2deg($distance); $distance = $distance * 60 * 1.1515; $distance = round($distance, 4); return $distance; } $slc_distance = distance_slc($lat_1, $lon_1, $lat_2, $lon_2); ?-->

PHP – Haversine Formula

<!--?php $earth_radius = 3960.00; # in miles $lat_1 = "47.117828"; $lon_1 = "-88.545625"; $lat_2 = "47.122223"; $lon_2 = "-88.568781"; $delta_lat = $lat_2 - $lat_1 ; $delta_lon = $lon_2 - $lon_1 ; function distance_haversine($lat1, $lon1, $lat2, $lon2) { global $earth_radius; global $delta_lat; global $delta_lon; $alpha = $delta_lat/2; $beta = $delta_lon/2; $a = sin(deg2rad($alpha)) * sin(deg2rad($alpha)) + cos(deg2rad($lat1)) * cos(deg2rad($lat2)) * sin(deg2rad($beta)) * sin(deg2rad($beta)) ; $c = asin(min(1, sqrt($a))); $distance = 2*$earth_radius * $c; $distance = round($distance, 4); return $distance; } $hav_distance = distance_haversine($lat_1, $lon_1, $lat_2, $lon_2); ?-->

Spherical Law of Cosines vs Haversine Formula

The table below lists the values of distances obtained using Spherical Law of Cosines and Haversine Formula, to many places from my favorite place on Earth.

 ##  Latitude1  Longitude1    Latitude2   Longitude2   SPHERICAL    HAVERSINE
 01  47.117828  -88.545625    47.122223   -88.568781      1.1302       1.1306
 02  47.117828  -88.545625    45.846967   -84.730038    201.6122     201.6853
 03  47.117828  -88.545625    42.487222   -83.083163    416.9887     417.1398
 04  47.117828  -88.545625    40.396333   -74.136083    853.9225     854.2319
 05  47.117828  -88.545625    14.230686    74.812998   8074.0414    8076.9674

Since most of our everyday measurements don’t need precision beyond 1st (or 2nd) decimal place, both these approaches are in excellent agreement with each other. As such, for my proposed API that will pull up locations nearest to a given pair of GPS coordinates from the local database to geotag my photos, I hope to use the Spherical Law of Cosines.

I understand that there could be other [more reliable] methods to precisely estimate the distance between two locations using their GPS coordinates and I [as well as other interested readers] would appreciate any such information.

22 Replies to “PHP – Calculating Distance Between Two Locations Given Their GPS Coordinates”

  1. Good one! Few days back,i wanted to build a web based system given GPS coordinates, but i wanted driving distance and not just the distance. I came up with two options:
    1. Using GoogleMaps API i can build a mashup. This is useful as i will have access to google map features. I wrote a simple javascript (thanks to Google and its API). No big deal infact 🙂
    2. Since all i needed was only driving distance, I used the following perl code:


    use LWP::Simple;
    use XML::Simple; # I’m just a simple guy.
    use Data::Dumper;

    $a = “16801”;
    $b = “61801”;

    my $f = get “{$a}+to+{$b}&output=kml”;
    $f =~ m/Distance: ([0-9.-]+)/;
    $distance = $1;
    print $distance;

  2. Adding to my last note: you can use LAT,LON in place of ZIP codes and you will still get driving distances

    Sorry for a small typo: $f =~ m/Distance: ([0-9.,-]+)/;

  3. Great formulars…

    I have even converted into a km. This formula is in miles as I think.

    However I have tried to insert this into a query which I run to get sights out of the DB and becase one function can be called once I am having a problems… Any solutions??


  4. Nice little forumla and explaination.
    In your SLC formular on line 19: $distance = $distance * 60 * 1.1515
    What is the “1.1515” supposed to represent?

    1. There are 60 nautical miles per degree separation of longitudes. And one nautical mile is same as 1.1515 statute (regular) mile. Since we measure distances in regular miles,

      $distance = $distance * 60 * 1.1515

      makes its appearance.

  5. Need some help answering the following question I read on…..I am ok with math, but not this much, lol Backround: There was a recent cloud formation over Russia that looked like a hole in the sky, and multiple major earthquakes on the Pacific plate. Here is the unanswered question from the site: Is there a correlation between ‘holes in the sky” and quakes? Too early to tell, that’s for sure. Let me see, here: Moscow is 55.75 North by 37.6 East and the quake epicenters can be found from USGS here, so with a little high powered math, if this hole in the sky stuff is real, we should be able to work the problem backwards (as with a bullet hole through a window and into an interior wall, right?) and figure out who is shooting at us.

    It would be interesting to see where this all points to.


  6. #!/usr/bin/perl
    use LWP::Simple;
    use XML::Simple; # I.m just a simple guy.
    use Data::Dumper;

    my $f = get “,-73.628195+to+45.462343,-73.628150&output=kml”;
    $f =~ m/Distance: ([0-9.,-]+)(.+m )/;
    $distance = $1;
    @v2 = split(/;/,$2);
    print $distance;
    print $um;

    print “\n\n”;

    my $f = get “,-73.628195+to+45.489785,-73.650043&output=kml”;
    $f =~ m/Distance: ([0-9.,-]+)(.+m )/;
    $distance = $1;
    @v2 = split(/;/,$2);
    print $distance;
    print $um;
    print “\n\n”;

  7. Great, now how do we take into account elevation … let’s say we one GPS coord is at the top of Everest, and the second GPS coord is one mile away by projection, but at the base of Everest… do we simply multiple by cos of the angle created by the two altitudes?

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.