Calculating distance between 2 points on Earth using the Haversine formula

20

Posted by Jaime | Posted in Coding | Posted on 17-02-2008

Tags: , ,

A while ago, I worked on a project where I had to implement code in C/C++ that could calculate the distance between 2 locations within the United States. It would be nice if the Earth was flat, since this would make my job a little easier, but the Earth is round and I needed to find a way that I could calculate distance. From the research that I have done, I found that plenty of people have worked on this problem before and one popular article is sited over and over again is this: R. W. Sinnott, “Virtues of the Haversine”, Sky and Telescope 68 (2), 159 (1984).

The Math used in all of the articles I read refer to the Haversine formula. Every example I found of this formula were either written in plain Math or shown in JavaScript. Although it was nice finding examples of this in JavaScript, there were a couple of problems with this:

  1. C datatypes aren’t the same as JavaScript datatypes and
  2. All examples used Latitude and Longitude, but I had to use Zip codes

So I had some more hurdles to overcome. They didn’t seem like big ones, but ones that I had to overcome nonetheless.

Now the funny thing about this project was that this was one of those endeavors where a spec wasn’t available before or during development. As a consequence, I would have designed the code differently than what I had have delivered. Not that there is anything wrong with what I ended up with, it’s just that I like neat code and what I wrote for this article is what I would have done if I had known all of the variables.

Example Code

What I am going to show you is how to calculate distance between two points within the United States, using the Haversine formula, implemented in C. The code will show you how to do that using Latitude and Longitude coordinates. The example application that is posted will show you how to use Zip codes in place of Latitude and Longitude.

The example code is demonstrated in an application written for Mac OS X. Win32 developers shouldn’t have to worry, since the Math code is done in a separate C file. Objective-C is used only for the UI file so you won’t have to read my Objective-C code to understand how to implement this formula in your own application. For the zip code implementation, I use the standard template library to contain the information I need in map so that I can extract the Latitude and Longitude for a particular postal code.

Creating the example code required these steps:

  1. Creating a Cocoa Application
  2. Creating the Haversine C code and
  3. Creating code that can handle Zip codes and finally
  4. Creating a routine to check the validity of the output

Creating a Cocoa Application

This is pretty simple step that most Apple developers are familiar with. I’m using Xcode 3.0, however, I set the project settings to be compatible with Xcode 2.4. Once the Cocoa application was created, I needed to create an interface that would take my Longitude and Latitude parameters and display the resulting output.

I also created a section that would take in Zip codes and get the result that way too! The logic is pretty simple and contained within the Objective-C code.

Creating the Haversine C code

This part was a little tough to work with. Every example I saw on the web used JavaScript. Here is the Math and JavaScript code I found at http://www.movable-type.co.uk/scripts/latlong.html :

Haversine Formula:
R = earth’s radius (mean radius = 6,371km)
Δlat = lat2− lat1
Δlong = long2− long1

a = sin²(Δlat/2) + cos(lat1) * cos(lat2) * sin²(Δlong/2)
c = 2 * atan2(√a, √(1−a))

d = R * c

(Note that angles need to be in radians to pass to trig functions).

JavaScript Implementation:

var R = 6371; // km
var dLat = (lat2-lat1).toRad();
var dLon = (lon2-lon1).toRad();
var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
 
Math.cos(lat1.toRad()) * Math.cos(lat2.toRad()) * Math.sin(dLon/2) * Math.sin(dLon/2);
 
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
var d = R * c;

Now for the C code!

To convert this code, there were three things I had to do to implement this as C code: A. Include the <math.h> file into my C file so that I get pi. B. I had to convert our Latitude and Longitude differences into radians, by multiplying the result by ( pi / 180 ) and C. I used the pow function to get the squared value instead of doing multiplying my value twice. Here is what I got:

double CalculateDistance( double nLat1, double nLon1, double nLat2, double nLon2 )
{
    int nRadius = 6371; // Earth's radius in Kilometers
    // Get the difference between our two points
    // then convert the difference into radians
    double nDLat = (nLat2 - nLat1) * (M_PI/180);
    double nDLon = (nLon2 - nLon1) * (M_PI/180);
    double nA = pow ( sin(nDLat/2), 2 ) + cos(nLat1) * cos(nLat2) * pow ( sin(nDLon/2), 2 );
 
    double nC = 2 * atan2( sqrt(nA), sqrt( 1 - nA ));
    double nD = nRadius * nC;
 
    return nD; // Return our calculated distance
}

Check out the application and you will see that it does NOT do the conversion properly. Here is where I had to do some digging around to see what was happening. The code on the website said one thing but the JavaScript code implemented on the website said another thing! I found that there were two lines omitted from the webpage: lat1 = lat1.toRad(), lat2 = lat2.toRad();. Of course, I had to add this to my code:

double CalculateDistance( double nLat1, double nLon1, double nLat2, double nLon2 )
{
    double nRadius = 6371; // Earth's radius in Kilometers
    // Get the difference between our two points
    // then convert the difference into radians
 
    double nDLat = ToRad(nLat2 - nLat1);
    double nDLon = ToRad(nLon2 - nLon1);
 
    // Here is the new line
    nLat1 =  ToRad(nLat1);
    nLat2 =  ToRad(nLat2);
 
    double nA = pow ( sin(nDLat/2), 2 ) + cos(nLat1) * cos(nLat2) * pow ( sin(nDLon/2), 2 );
 
    double nC = 2 * atan2( sqrt(nA), sqrt( 1 - nA ));
    double nD = nRadius * nC;
 
    return nD; // Return our calculated distance
}

Anyone who is looking the JavaScript code will wonder where to get the toRad function. This was actually declared within  the webpage, and I implemented my own C version called ToRad.

Now let’s move onto incorporating Zip codes into our example!

Creating code that can handle Zip codes

Now, to be able to use Zip codes in our example means that we have to be able to map our zip code to a Latitude/Longitude pair. Fortunately for us, there is a complete listing of postal codes at http://sourceforge.net/projects/zips/ , complete with town, state and postal code. All we have to do, is to read the values into our application and then create a routine for getting our Latitude/Longitude with a given Zip code.

For those who don’t know, Zip is actually an acronym for Zoning Improvement Plan. If you do a little research, you will find that there is an interesting history behind this 5 digit number!

In order for us to use our text file, there are a couple of things we have to do:

  1. Read the text file contents into memory.
  2. Parse the file contents into a container variable.
  3. Create a function that allows us to easy obtain the Latitude and Longitude for a given Zip code.

Since the above listed steps can be accomplished using cross platform friendly POSIX code, I will implement the C code in a C++ file to keep things nice and neat. The reason for putting my code in a C++ file is because I wanted to create a C++ struct that would hold my Latitude Longitude pair. Using a C++ struct allows me to create private and public struct members and even a constructor.

Create a routine to check the validity of the output

This part was pretty easy. Since the JavaScript implementation works, I created a simpler version of the JavaScript code and posted that page at the end of this article. Using that code, and Google, I can see that the application is outputting the right distance.

I hoped you enjoyed reading this article as much as I did writing it ;) Check out the sample code and happy coding!

Haversine application for Mac OS X

JavaScript Implementation of Haversine Formula

C Source Code and Xcode Project

Comments (20)

Thanks! I found the same example of javascript and was going to convert it to something C based, you’ve saved me some hassle. Excellent write-up as well. =)

Ugh.

So I’m doing some coding for the iPhone. I got everything written up nice and fun for the iPhone… the only problem is that apple’s already done this for me with this function:

getDistanceFrom:
Returns the distance (in meters) from the receiver’s coordinate to the coordinate of the specified location.

Oh well, it was an interesting process. Thanks for the education. =)

I’m glad I was able to help! It’s also very nice of Apple to give you a method that provides the functionality you needed!

When I first started working on .NET many moons ago, the first version of .NET did not have FTP support and I needed to write an application with .NET that sent files up to a FTP server. Who would figure that a technology that is Internet centric would not have basic FTP support.

It’s nice to see that Apple didn’t follow the same route as Microsoft.

-Jaime

[...] most popular post is actually the Haversine formula [...]

Yes, nice work, the explanation is pretty simple and powerful. Thanx

Glad I was able to help!

Sorry for the dumbness, but the distance returned in the method is in what ? Kilometers ? Tks for the help ..

The distance is returned to you in Kilometers.

Just for your information,

Earth’s radius in Kilometers
6372.797

for better precision.

Thanks for taking the time to publish this.
Paul K.

You’re welcome! It took me a while to put this together but it was definitely worth it.

“2.atan2(?a, ?(1?a))” can be written simpler (for all arguments) as “2.asin(?a)”

As Bob points out, using the actual haversine formula makes things a little simpler: you can replace

double nA = pow ( sin(nDLat/2), 2 ) + cos(nLat1) * cos(nLat2) * pow ( sin(nDLon/2), 2 );
double nC = 2 * atan2( sqrt(nA), sqrt( 1 – nA ));

with the slightly simpler

double nA = haversin(nDLat) + cos(nLat1) * cos(nLat2) * haversin(nDLon);
double nC = 2 * asin( sqrt(nA) );

I have no idea what Kiichi is talking about.
The Wikipedia: Earth radius article points out that the radius various from 6,357 km to 6,378 km, and suggests that if you must choose one particular value, “6371 km” is a good mean value.

Do you know if with little improvement you can calculate distance by sea only using lat/long for A/B points ?

[...] found a helpful post from Jaime Rios explaining the Haversine Formula. MySql provides the necessary functions to translate Jaime’s work into SQL. The Haversine [...]

Hi Boris,
If you look at the function CalculateDistance, you see that it takes a latitude and longitude. @boris

Thanks for the function, it worked out great for me with some small changes to the parameters but the core function was what I needed.

This was a very interesting, and rewarding, article for me to write, so I’m glad that it was able to help.

Do you happen to have a code using the Haversine (or know of one) That actually asks the user in C++ to input the longitude and latitude themselves to find the distance in miles?

That is a pretty simple console application to make. Are you using that for a homework assignment or something?

Write a comment