On Thu, 14 Jan 2010 10:27:17 +0100 Karl Ove Hufthammer <k...@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 _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
