# 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

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

1 <!--?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”

1. Vikas Argod says:

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; 2. Vikas Argod says: 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??

Regards,

1. I can try to take a look at it – can you please send me [parts of] the query where you are using this ?

Regards,

3. Good to hear [but sorry for delayed response from my end though] 🙂

4. 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

5. Robert says:

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.

6. Krzyzmo says:

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

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

8. chamara says:

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

9. 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!

10. #!/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”;

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

12. Joe says:

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?

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