Distance and Bearing Between Latitude/Longitude Points

September 16th, 2010

A recent project required me to figure out how to calculate the distance and bearing between two latitude/longitude points. It was a bit of a challenge, but I finally got it working. Here’s how it’s done.

It’s important to keep in mind that these formulas work in radians, not degrees. So if you are implementing these equations with what I consider “normal” lat/long pairs in degrees, such as (32.853135, -96.971728), you must convert the values to radians before performing the calculation. Note the conversions in my example implementations.

Distance

There are two equations for calculating distance between two latitude/longitude points. The first is called the Haversine Formula:

R = Earth’s radius (6,371km, or 3,959mi)
Δ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

A slightly easier to implement version of the calculation is called the Vincenty Formula:

d = acos(sin(lat1).sin(lat2)+cos(lat1).cos(lat2).cos(long2−long1)).R

Expressed as a Microsoft SQL Server 2008 function, the Vincenty Formula looks like this:

CREATE FUNCTION dbo.Distance
(
@inlat1 decimal(9,6),
@inlong1 decimal(9,6),
@inlat2 decimal(9,6),
@inlong2 decimal(9,6)
)
RETURNS decimal(5,1)
AS
BEGIN
declare @return decimal(5,1)

declare @lat1 float
declare @long1 float
declare @lat2 float
declare @long2 float

declare @x float
declare @y float

select @lat1=radians(@inlat1)
select @long1=radians(@inlong1)
select @lat2=radians(@inlat2)
select @long2=radians(@inlong2)

select @return=Acos(sin(@lat1)*sin(@lat2)+cos(@lat1)*cos(@lat2)*cos(@long2-@long1)) * 3959

RETURN @return
END

Bearing

To calculate the intial bearing, or direction, between the two points, the formula is:

θ = atan2(sin(Δlong).cos(lat2), cos(lat1).sin(lat2) − sin(lat1).cos(lat2).cos(Δlong))

I implemented the bearing function in Microsoft SQL Server 2008, where there is no atan2 function, like this:

CREATE FUNCTION dbo.Bearing
(
@inlat1 decimal(9,6),
@inlong1 decimal(9,6),
@inlat2 decimal(9,6),
@inlong2 decimal(9,6)
)
RETURNS decimal(3,0)
AS
BEGIN
declare @return decimal(3,0)

declare @lat1 float
declare @long1 float
declare @lat2 float
declare @long2 float

declare @x float
declare @y float

select @lat1=radians(@inlat1)
select @long1=radians(@inlong1)
select @lat2=radians(@inlat2)
select @long2=radians(@inlong2)

select @y=sin(@long2-@long1)*cos(@lat2)
select @x=cos(@lat1)*sin(@lat2)-sin(@lat1)*cos(@lat2)*cos(@long2-@long1)

if(@y>0) begin
if(@x>0) begin
select @return=degrees(atan(@y/@x))
end
if(@x<0) begin
select @return=180-degrees(atan(-@y/@x))
end
if(@x=0) begin
select @return=90
end
end
if(@y<0) begin
if(@x>0) begin
select @return=-degrees(atan(-@y/@x))
end
if(@x<0) begin
select @return=degrees(atan(@y/@x))-180
end
if(@x=0) begin
select @return=270
end
end
if(@y=0) begin
if(@x>0) begin
select @return=0
end
if(@x<0) begin
select @return=180
end
if(@x=0) begin
select @return=0
end
end

select @return=@return+180

RETURN @return
END

This returns a bearing in positive degrees with 0 equal to North, rounded to 0 decimal places.

Destination

The next logical question is, how do I find a destination point given a start point and initial bearing. Here's the equation for that:

lat2 = asin(sin(lat1)*cos(d/R) + cos(lat1)*sin(d/R)*cos(θ))
lon2 = lon1 + atan2(sin(θ)*sin(d/R)*cos(lat1), cos(d/R)−sin(lat1)*sin(lat2))

I leave the implementation of that one to you!