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 and are the Cartesian coordinates of two points, the distance between them, , would be:

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 is the radius of our beautiful Earth, the elevation angle (almost similar to *latitude* – almost, yeah; explanation comes soon) and the azimuthal angle (similar to longitude), then we can write the points and as and . The transformation from Cartesian to Spherical Polar coordinates is not that difficult and may be written as

If the straight lines drawn through these points, and , intersect each other at the center of the Earth and subtend an angle , then by extending the Law of Cosines [from plane trigonometry] to the spherical case, it can be shown that

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

It is generally assumed that elevation angle, , can take values in the range and azimuthal angle, , can take values in the range . However, the accepted geographic convention says that latitude (elevation angle) can range between and that longitude (azimuth angle) can range between . The expression we have above for is still valid, except that we need to replace with and by . Once those replacements are made, the result becomes

where indeed represents the latitude and 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 is expressed in degrees (not radians), then the distance between points and will be

This looks nice, elegant and easily computable but has an inherent problem: since is not so well behaved, the resulting distance could be erroneous if the longitudinal separation between and 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,

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

The in last step is a result of the fact that 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 and are antipodal – this does provide a pretty accurate estimate for distance.

### PHP – Spherical Law of Cosines

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 | <?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 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 | <?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.

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;

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

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,

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

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?

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.

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

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,

Good to hear [but sorry for delayed response from my end though] ðŸ™‚

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

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.

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

thank you

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!

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

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?

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?