Converting Latitude and Longitude into Distance

Solving for Great Circle distance

Calculating distance when using flat map data using a rectangular grid coordinates system is a simple matter of using Pythagoras to solve for the distance:

`    'Compute the components and the distance...`
`    myDeltaX = Abs(myOriginX - myDestX)`
`    myDeltaY = Abs(myOriginY - myDestY)`
`    myDistance = Sqr((myDeltaX ^ 2) + (myDeltaY ^ 2))`
` `

It different with latitude and longitude. For very short distances (less than about 2km) you can get away with using pythagoras to calculate the distance between latitude and longitude points. The error is small - less in fact than the error that likely is in your data already.

The following functions will solve for the great circle distance between two points where the coordinates are specified in decimal degrees of latitude and longitude, assuming no change in elevation. Longitudes in the western hemisphere must be specified as negative values, in the eastern hemisphere, they are positive. Use negative latitudes in the southern hemisphere, use positive latitudes the north. Coordinates must be entered in decimal degrees of longitude and latitude. Longitude for the western hemisphere and latitude for the southern hemisphere are expressed as negative values.

The function "LatLonDistance" has the following parameters:

1. Latitude of point 1 - Double
2. Longitude of point 1 - Double
3. Latitude of point 2 - Double
4. Longitude of point 2 - Double
5. Units - string

Units must be one of the following:

1. miles - "mi"
2. feet - "ft"
3. yards - "yd"
4. kilometers - "km"
5. meters - "m"

The function returns the distance (in the specified units) as a Double. If invalid or no units are specified, the function returns -1.

by Chris North, January 1997

Visual Basic 3.0 and 4.0 Source Code:

 Remember: The underscore line continuation character ("_") will only work in Visual Basic 4.0. You will have to remove it and have those lines of code here that end with the ("_") character all on one line for Visual Basic 3.0
`Function LatLonDistance (ByVal pDb_Lat1 As Double, _`
`                         ByVal pDb_Lon1 As Double, _`
`                         ByVal pDb_Lat2 As Double, _`
`                         ByVal pDb_Lon2 As Double, _`
`                         ByVal pSt_Units As String) As Double`
`Dim lDb_R As Long 'Radius of Earth = 6367 km = 3956 mi`
`Dim lDb_deltaLat As Double`
`Dim lDb_deltaLon As Double`
`Dim lDb_a As Double`
`Dim lDb_c As Double`
` `
`    'Set the radius of the earth in the desired units...`
`    Select Case UCase(pSt_Units)`
`    Case "MI" ' Miles`
`        lDb_R = 3956`
`    Case "FT" ' Feet`
`        lDb_R = 20887680`
`    Case "YD" ' Yards`
`        lDb_R = 6962560`
`    Case "KM" ' Kilometres`
`        lDb_R = 6367`
`    Case "M"  ' Metres`
`        lDb_R = 6367000`
`    Case Else ' non-supported units`
`        LatLonDistance = -1`
`        Exit Function`
`    End Select`
` `
`    'Calculate the Deltas...`
`    lDb_deltaLon = AsRadians(pDb_Lon2) - AsRadians(pDb_Lon1)`
`    lDb_deltaLat = AsRadians(pDb_Lat2) - AsRadians(pDb_Lat1)`
` `
`    'Intermediate values...`
`    lDb_a = Sin2(lDb_deltaLat / 2) + _`
`        Cos(AsRadians(pDb_Lat1)) * _`
`        Cos(AsRadians(pDb_Lat2)) * _`
`        Sin2(lDb_deltaLon / 2)`
`    `
`    'Intermediate result c is the great circle distance in radians...`
`    lDb_c = 2 * Arcsin(GetMin(1, Sqr(lDb_a)))`
` `
`    'Multiply the radians by the radius to get the distance in specified units...`
`    LatLonDistance = lDb_R * lDb_c`
` `
`End Function`
` `
`Private Function Arcsin (ByVal X As Double) As Double`
`    'Arcsin(X) = Atn(X / Sqr(-X * X + 1)) [from MS Help - VB4, 1995]`
`    Arcsin = Atn(X / Sqr(-X * X + 1))`
`End Function`
` `
`Private Function AsRadians (ByVal pDb_Degrees As Double) As Double`
`Const vbPi = 3.14159265358979`
`    'To convert decimal degrees to radians, multiply`
`    'the number of degrees by pi/180 = 0.017453293 radians/degree`
`    AsRadians = pDb_Degrees * (vbPi / 180)`
`End Function`
` `
`Private Function GetMin (ByVal X As Double, ByVal Y As Double) As Double`
`'The min() function protects against possible roundoff errors that`
`'could sabotage computation of the arcsine if the two points are`
`'very nearly antipodal (that is, on opposide sides of the Earth)`
`'                                       - RGC, 1996`
`    If X <= Y Then`
`        GetMin = X`
`    Else`
`        GetMin = Y`
`    End If`
`End Function`
` `
`Private Function Sin2 (ByVal X As Double) As Double`
`    'sin^2(X) = (1 - cos(2X))/2 [from Greer and Hancox, 1991]`
`    Sin2 = (1 - Cos(2 * X)) / 2`
`End Function`

The best way to calculate the distance between 2 points on a globes surface.

The Earth is round, but big, so we can consider it flat for short  distances. But, even though the circumference of the Earth is about  25,000 miles (40,000 kilometers), flat-Earth formulas for calculating  the distance between two points start showing noticeable errors when  the distance is more than about 12 miles (20 kilometers). Of course,  how much error is "noticeable" depends on how you are going to use the  result.

Cartesian coordinates express distances in two different directions,  such as north-south for one direction and east-west for the other.  The straight line distance between two points can then be thought of  as the long side of a right triangle with one of the short sides being  the north-south distance between the points and the other being the  east-west distance. (A right triangle is one that has a square  corner.) The usual formula for computing the length of the long side  of a right triangle is the Pythagorean Theorem. Using this formula  from geometry requires knowing about square roots.

Near the North Pole and near the South Pole, the longitude lines,  which go north-south and are called the meridians, approach each other  noticeably - in fact, they meet at the pole. The latitude lines, which  go east-west, are circles around the pole. Treating differences in  locations along these directions as if they were the sides of a right  triangle leads to errors in the computation of distance. Very close  to the pole, the answer could be VERY wrong - but a different  flat-Earth approximation, obtained from plane trigonometry, can be  used for short distances: the Polar Coordinate Flat-Earth Formula.  Using this formula - and the others mentioned below - requires knowing  something about trigonometry.

When you study spherical trigonometry, you learn a bunch of formulas.  One of them is called the Law of Cosines for Spherical Trigonometry.  It is a perfectly fine formula when it is used for the right purposes.  Solving for short distances on the surface of the Earth is not one of  those purposes. The problem is as follows: Suppose you have a right  triangle with a very small angle. The ratio between the short side and the long side is very close to 1.0. The formula  computes that ratio first, then requires the computer to find the angle  that has that ratio. In principle, the computer can do so but ordinary computers  approximate all numbers to some number of what are called "significant  digits". With 7 or 8 significant digits, the computer cannot  distinguish between the ratios for angles smaller than about a minute  of arc (a minute is 1/60 of a degree). Since the angle being computed  has its apex at the center of the Earth, a minute of arc corresponds to  about a mile on the surface.

Since the formula is mathematically correct, it can be manipulated  into other forms. The Haversine Formula is one result of such  manipulations. It has a similar problem, but it is "poorly  conditioned" when the two points are all the way around the Earth  from each other, rather than when they are close to each other. The  discussion below gives a second version of the haversine formula that  is easier to program on some computers.

It also discusses the fact that the Earth is not quite a sphere, and  gives several references for further reading. There is discussion of  where to look if even the ellipsoid approximation is too coarse. It  also has a pointer to one source of navigation formulas on the  Internet.

~~~~~~~~~~~~~~~~~~~~~~

From: Bob Chamberlain <rgc@jpl.nasa.gov>

The distances considered here are along the surface of the Earth,  and ignore the effect of differences in elevation.

If the distance is less than about 20 km (12 mi) and the locations of  the two points in Cartesian coordinates are X1,Y1 and X2,Y2 then the Pythagorean Theorem:

d = sqrt((X2 - X1)^2 + (Y2 - Y1)^2)

will work well for latitudes less than 70 degrees and have an error of less than 30 meters (100 ft) , and less than 20 meters ( 66 ft) for latitudes less than 50 degrees, and less than 9 meters ( 30 ft) for latitudes less than 30 degrees  (These error statements reflect both the convergence of the meridians  and the curvature of the parallels. The error is non-linear with  distance; shorter distances will have better percentage errors.)

The flat-Earth distance d will be expressed in the same units as the  coordinates.

Presuming a spherical Earth with radius R, and  the locations of the two points in spherical coordinates (longitude  and latitude) are lon1,lat1 and lon2,lat2 then theHaversine Formula:

dlon = lon2 - lon1
dlat = lat2 - lat1
a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
c = 2 * arcsin(min(1,sqrt(a)))
d = R * c

will give mathematically and computationally exact results. The  intermediate result c is the great circle distance in radians. The  great circle distance d will be in the same units as R.

When the two points are antipodal (on opposite sides of the Earth),  the Haversine Formula is ill-conditioned, but the error, perhaps  as large as 2 km (1 mi), is in the context of a distance near  20,000 km (12,000 mi). Further, there is a possibility that round-off  errors might cause the value of sqrt(a) to exceed 1.0, which would  cause the inverse sine to crash without the bulletproofing provided by  the min() function.

Most computers require the arguments of trigonometric functions to be  expressed in radians. To convert lon1,lat1 and lon2,lat2 from degrees,  minutes, and seconds to radians, first convert them to decimal  degrees. To convert decimal degrees to radians, multiply the number  of degrees by pi/180 = 0.017453293 radians/degree.

Inverse trigonometric functions return results expressed in radians.  To express c in decimal degrees, multiply the number of radians by  180/pi = 57.295780 degrees/radian. (But be sure to multiply the number  of RADIANS by R to get d.)

The Haversine Formula can be expressed in terms of a two-argument  inverse tangent function, atan2(y,x), instead of an inverse sine  as follows:

dlon = lon2 - lon1
dlat = lat2 - lat1
a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
c = 2 * atan2( sqrt(a), sqrt(1-a) )
d = R * c

In this context, a one-argument inverse tangent function would also  work: replace "atan2( sqrt(a), sqrt(1-a) )" by "arctan( sqrt(a) /  sqrt(1-a) )", but insert a test to ensure you will not divide by  zero if a = 1. (In the case a = 1, c = pi radians = 180 degrees,  and d is halfway around the Earth = 3.14159265... * R.)

The problem of determining the great circle distance on a sphere has  been around for hundreds of years, as have both the Law of Cosines  solution and the Haversine Formula.  Sinnott gets the credit here because he was quoted by Snyder.

The Pythagorean flat-Earth approximation assumes that meridians are  parallel, that the parallels of latitude are negligibly different from  great circles, and that great circles are negligibly different from  straight lines. Close to the poles, the parallels of latitude are not  only shorter than great circles, but indispensably curved. Taking this  into account leads to the use of polar coordinates and the planar law  of cosines for computing short distances near the poles the

Polar Coordinate Flat-Earth Formula:

a = pi/2 - lat1
b = pi/2 - lat2
c = sqrt( a^2 + b^2 - 2 * a * b * cos(lon2 - lon1) )
d = R * c

is computationally only a little more expensive than the Pythagorean  Theorem and will give smaller maximum errors for higher latitudes and  greater distances. The maximum errors, which depend upon azimuth in  addition to separation distance, are equal at 80 degrees latitude when  the separation is 33 km (20 mi), 82 degrees at 18 km (11 mi),  84 degrees at 9 km (5.4 mi). But even at 88 degrees the polar error  can be as large as 20 meters (66 ft) when the distance between the  points is 20 km (12 mi).

The latitudes lat1 and lat2 must be expressed in radians;  pi/2 = 1.5707963. Again, the intermediate result c is the distance in  radians and the distance d is in the same units as R.

An UNRELIABLE way to calculate distance on a spherical Earth is the

Law of Cosines for Spherical Trigonometry:

a = sin(lat1) * sin(lat2)
b = cos(lat1) * cos(lat2) * cos(lon2 - lon1)
c = arccos(a + b)
d = R * c

Although this formula is mathematically exact, it is unreliable for  small distances because the inverse cosine is ill-conditioned. Sinnott  (in the article cited above) offers the following table to illustrate  the point:

cos (5 degrees) = 0.996194698
cos (1 degree) = 0.999847695
cos (1 minute) = 0.9999999577
cos (1 second) = 0.9999999999882
cos (0.05 sec) = 0.999999999999971

A computer carrying seven significant figures cannot distinguish the  cosines of any distances smaller than about one minute of arc.  The function min(1,(a + b)) could replace (a + b) as the argument for  the inverse cosine to guard against possible round-off errors, but  doing so would be to "polish a cannonball".

~~~~~~~~~~~~~~~~~~~~~~

What value should I use for the radius of the Earth, R?

The historical definition of a "nautical mile" is "one minute of arc  of a great circle of the earth". Since the earth is not a perfect  sphere, that definition is ambiguous. However, the internationally  accepted (SI) value for the length of a nautical mile is (exactly, by  definition) 1.852 km or exactly 1.852/1.609344 international miles  (that is, approximately 1.15078 miles - either "international" or  "U.S. statute"). Thus, the implied "official" circumference is 360  degrees times 60 minutes/degree times 1.852 km/minute = 40003.2 km.  The implied radius is the circumference divided by 2 pi:

R = 6367 km = 3956 mi

~~~~~~~~~~~~~~~~~~~~~~

When is it NOT okay to assume the Earth is a sphere?

A quick test is: Compare the results produced by using the two  extreme values of the radius of curvature for the Earth:

minimum radius of curvature: 6336 km (3937 mi)
maximum radius of curvature: 6399 km (3976 mi)

in your application. If the results are different enough to cause you  to change your action (or your recommendation, or your interpretation  of the implication of the results, etc.), then assuming the Earth is  spherical is NOT okay.

The shape of the Earth is well approximated by an oblate spheroid. The  radius of curvature varies with direction and latitude. According to  formulas given on pages 24 and 25 of the book by Snyder,

The radius of curvature of an ellipsoidal Earth in the plane of the  meridian is given by

R' = a * (1 - e^2) / (1 - e^2 * (sin(lat))^2)^(3/2)

where a is the equatorial radius,  b is the polar radius, and  e is the eccentricity of the ellipsoid = sqrt(1 - b^2/a^2)
and the radius of curvature in a plane perpendicular to the meridian  and perpendicular to a plane tangent to the surface is given by

N = a/sqrt(1-e^2*(sin(lat)^2))

A Swedish book  suggests use of the geometric mean of these two radii of  curvature for all azimuths, as it produces errors of order of  magnitude 0.1% for distances within 500 km (300 mi) at 60 degrees  latitude. The formula for that average is no more complicated than  either of its components. That is, r = sqrt(R' * N) becomes

r = a * sqrt(1 - e^2) / (1 - e^2 * (sin(lat))^2)

Using these formulas with

a = 6378 km (3963 mi) Equatorial radius (surface to center distance)
b = 6357 km (3950 mi) Polar radius (surface to center distance)
e = 0.081082 Eccentricity

gives the following table of values for the Radii of Curvature:

latitude...........r...................R'..................N..........
00 degrees . 6357 km (3950 mi) . 6336 km (3937 mi) . 6378 km (3963 mi)
15 degrees . 6360 km (3952 mi) . 6340 km (3940 mi) . 6379 km (3964 mi)
30 degrees . 6367 km (3957 mi) . 6352 km (3947 mi) . 6383 km (3966 mi)
45 degrees . 6378 km (3963 mi) . 6367 km (3957 mi) . 6389 km (3970 mi)
60 degrees . 6388 km (3970 mi) . 6383 km (3966 mi) . 6394 km (3973 mi)
75 degrees . 6396 km (3974 mi) . 6395 km (3974 mi) . 6398 km (3975 mi)
90 degrees . 6399 km (3976 mi) . 6399 km (3976 mi) . 6399 km (3976 mi)

Note that the radius of curvature for an ellipsoid is not the same as  the distance from the surface of the ellipsoid to the center. In fact,  the radius of curvature increases as the radius decreases. Also, be  aware that a variety of ellipsoids with slightly different parameters  have been fit to the Earth; the preferred ellipsoid may depend on the  region in which you are most interested.

The  spherical earth computations will provide  underestimates of real world distances measured in the direction of  the equator (and especially for trans-equatorial links) and  overestimates for those measured in the direction of the poles (and  especially for trans-polar ones).

For most purposes, it is quite satisfactory to treat the Earth as a  sphere. If not, an ellipsoid can provide a better approximation.

Software for solving distance and azimuth problems on the ellipsoid can be obtained (as of 10/10/96) by anonymous ftp from several sources, two of which are listed below:

The URL of the National Geodetic Survey (of the National Oceanic and
Atmospheric Administration in the US Department of Commerce) is:

Another anonymous ftp source for ellipsoid software is the US
Geological Survey (of the US Department of the Interior), at:

~~~~~~~~~~~~~~~~~~~~~~

When is it NOT okay to assume the Earth is an ellipsoid

The shape the Earth would assume if it were all measured at mean sea level is called the geoid. The geoid varies no more than about a hundred meters above or below a well-fitting ellipsoid, a variation far less than the ellipsoid varies from the sphere. Terrain relief is reported relative to the geoid. (Paraphrased from p. 11 of the book by Snyder cited above.)

Distances on the surface of the geoid are not particularly meaningful. However, there are applications, such as long-term prediction of orbits of Earth satellites, that require better approximations than are provided by an ellipsoid.

~~~~~~~~~~~~~~~~~~~~~~

Where can I find formulas for great circle navigation problems,

like what course to follow and where will I be after following a known azimuth for a given distance?

See Ed Williams' navigation formulas at

which can also be accessed by clicking on 'FAQ Lists' at

~~~~~~~~~~~~~~~~~~~~~~

How can I calculate distance on an ellipsoid.

I found two approaches to this problem.

One approach involves computing a discretized approximation to a geodesic. This works on any ellipsoid. The first two software pointers follow this approach. As far as I can tell, all the other references follow the second approach.

The second approach is designed for the Earth, exploiting two special features: (1) it is an ellipsoid of revolution, and (2) it is nearly spherical. These properties allow one to construct a solution as a power series in the small "flattening" parameter. One truncates the series depending upon the accuracy required.

I was looking for a solution that works for any ellipsoid. Most of what follows is geodesy-centric, so I have only explored a handful of the suggestions listed below. Also --- as one person noted --- a web search with key words of "ellipsoid & distance & arc" will generate a list of hundreds of sites. (All dealing with the Earth, I might add).

What follows is heavily cut'n'pasted from the responses or from web pages. Don't take these as advice or endorsement from me -- these are not my words below!

Available Software on the net:

The code here computes a polygonal approximation to a geodesic and then refines it until the geodesic curvature vanishes at each vertex on the path.

Similar to the above, except that it actually works on any manifold --
you supply a routine to compute the metric at any point, etc.

Well-proven code, forward and inverse problem.

>From the (U.S.) National Geodetic Survey:

INVERSE/FORWARD3D (Version 1.0)
Comprises four programs - Inverse (Version 2.0) which computes the  geodetic azimuth and distance between two points, given their  geographic positions; Forward (Version 2.0) which computes the  geographic position of a point, given the geodetic azimuth and  distance from a point with known geographic position; and the  three-dimensional versions of these programs . INVERS3D (Version 1.0)  and FORWRD3D (Version 1.0), which include the height component.

This program will solve either the geodetic position or inverse  problems in any quadrant of the Earth using the parameters of one of  five commonly-used ellipsoids.

A product called NACNav (60USD/license) which can calculate the  shortest distance between two points anywhere on the earth surface and  the angle between the direction and the magnetic north (determined by  local magnetic declination).

========================
This answer was prepared by Robert G. Chamberlain of Caltech (JPL),  <rgc@jpl.nasa.gov>, and reviewed on the comp.infosystems.gis newsgroup  in Oct 1996. It was revised in November 1997 and February 1998.  Formatting was updated and the hyperlinks were checked in January  2000. Minor cosmetic changes were made and hyperlinks were updated  in February 2001.

Steven Michael Robbins' discussion of distance computation on an ellipsoid was added in May 1999.

The pre-geometry summary was added in January 2000 in response to a request from Karen Kast, who was seeking information that might be useful to a 6th grade student doing a science project.