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}


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

Related Posts