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'
 
) CHARACTER SET utf8 ;

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}$

While 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 Geonames.org 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 Geonames.org 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”

Vikas Argod says:

Tuesday, 2009-08-04 at 3:42 pm -0400

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:

#!/usr/bin/perl

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

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

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

Vikas Argod says:

Tuesday, 2009-08-04 at 3:44 pm -0400

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.,-]+)/;

Pingback: Gowtham :: Brain Dump v1.0 Â» PHP â€“ Calculating Distance Between …
Pingback: Gowtham :: Brain Dump v1.0 » MySQL – Finding Locations Nearest To A Given Pair Of GPS Coordinates
Goran says:

Thursday, 2009-08-20 at 10:30 am -0400

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??

Regards,

1. Gowtham says:
  
  Tuesday, 2009-09-08 at 7:49 am -0400
  
  I can try to take a look at it – can you please send me [parts of] the query where you are using this ?
2. Goran says:
  
  Tuesday, 2009-09-08 at 7:51 am -0400
  
  10x for the reply, but I’ve already manage it.
  
  take a look: http://travel2macedonia.com.mk/travel-to-macedonia/destinations/skopje/sights-to-visit
  
  Regards,
3. Gowtham says:
  
  Tuesday, 2009-09-08 at 7:53 am -0400
  
  Good to hear [but sorry for delayed response from my end though] :)
Debt Settlement Program says:

Wednesday, 2009-09-02 at 2:20 pm -0400

Apologize for my bad english, I over its a winsome drama of your writing. Sumptuously I be suffering with faced alot of difficulties in this condition but your article will definately eschew me in future. Say thank you You

Robert says:

Monday, 2009-09-07 at 9:16 pm -0400

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. Gowtham says:
  
  Tuesday, 2009-09-08 at 7:48 am -0400
  
  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.
Krzyzmo says:

Friday, 2009-10-09 at 6:06 pm -0400

Need some help answering the following question I read on urbansurvival.com…..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.

Krzyzmo

back care products says:

Tuesday, 2009-12-15 at 10:20 pm -0500

This is an spectacular entry. Thank you very much for the supreme post provided! I was looking for this entry for a long time, but I wasn’t able to find a reliable source.

chamara says:

Friday, 2010-02-12 at 4:22 am -0500

I want to calculate distance between two pints using AWK programming.If you know about it please send me amessege
thank you

Leann Musselwhite says:

Wednesday, 2010-03-10 at 5:56 pm -0500

Really great article – I was heading for a similar article which I will probably still write, but from a slightly different angle. Thanks for sharing this with your readers…Obviously a lot of others appreciate it too!

cristian l says:

Monday, 2010-05-24 at 2:16 am -0400

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

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

print “\n\n”;

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

Pingback: Ordnance Survey 1:50 000 Gazetteer import and reuse « Central Unit blogs
Pingback: Ordnance Survey 1:50 000 Gazetteer import and reuse « Portable Antiquities Scheme blogs
Tom Delagado says:

Wednesday, 2010-07-14 at 10:55 am -0400

Hello just thought i would tell you something.. This is twice now i’ve landed on your blog in the last 3 days searching for totally unrelated points. Spooky or what?

Pingback: Dichtstbijzijnde adressen zoeken via geocoÃ¶rdinaten in PHP « RRWD studie webtoepassingen
Pingback: links for 2011-02-21 « Kaigani’s Arbor Vitae
Joe says:

Monday, 2011-09-26 at 11:19 am -0400

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?

Underlying Mathematics

PHP – Spherical Law of Cosines

PHP – Haversine Formula

Spherical Law of Cosines vs Haversine Formula

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

Leave a Reply Cancel reply