[R-sig-Geo] Best way to calculate distance between geographical points

[Karl Ove Hufthammer]

On Thu, 14 Jan 2010 10:27:17 +0100 Karl Ove Hufthammer <karl at huftis.org> 
wrote:
> I have found three functions for calculating the distance between 
> geographical points, 'geo.dist' in the 'oce' package, 'rdist.earth' in 
> the 'fields' package, and 'spDistsN1' in the 'sp' package. They all give 
> slightly difference answers, so my question is: which of them is best, 
> i.e., gives the most accurate answer?

One simple test would be to calculate the circumference of the earth at 
the equator and the meridian. But this test gives conflicting answers.

According to Wikipedia, "http://en.wikipedia.org/wiki/Earth" the 
equatorial circumference is 40075.02 km and the meridional circumference 
is 40007.86 km.

When calculating half the equatorial circumference, 'spDistsN1' is best, 
and gives the *exact* answer, while 'rdist.earth' is close (less than 
one km difference), and 'geo.dist' is far from the correct answer 
(133.92 km difference). So it would initially seem that 'spDistsN1' is 
the one to use, and one should avoid 'geo.dist'.

But when calculating the meridional distance, 'geo.dist' gives the exact 
answer, while both 'rdist.earth' and 'spDistsN1' are both off by about 
34 km (but in difference directions).

So which to use is still an open question ...

Here is the code I used:

library(oce)
library(fields)
library(sp)

# Equatorial
lt1=0
lt2=0
lg1=0
lg2=180

a=matrix(c(lg1,lt1),nrow=1)
b=matrix(c(lg2,lt2),nrow=1)

# Calculated distance
geod.dist(lt1,lg1,lt2,lg2)
rdist.earth(a,b, miles=FALSE)
spDistsN1(a,b, longlat=TRUE)

# Correct distance (according to Wikipedia)
40075.02/2


# Meridional
lt1=-90
lt2=90
lg1=0
lg2=0

a=matrix(c(lg1,lt1),nrow=1)
b=matrix(c(lg2,lt2),nrow=1)

# Calculated distance
geod.dist(lt1,lg1,lt2,lg2)
rdist.earth(a,b, miles=FALSE)
spDistsN1(a,b, longlat=TRUE)

# Correct distance (according to Wikipedia)
40007.86/2 

-- 
Karl Ove Hufthammer