1

I am using the function distVincentyEllipsoid() from the geosphere package to calculate distances between the world's ports. This function takes as input the longitude and latitude of two places and gives the distance in meters as output. It appears to work well generally, but in some cases it's been giving me NAs for no discernible reason. I was wondering whether anyone was aware of a bug (there is no documentation online I could find that addresses this issue), or could perhaps point me to what I'm doing wrong. What's strange is that it appears to affect only specific pairs – eg. f(a,b) and f(b,a) will give NA, but f(a,c) and f(b,c) will not.

Eg. Whangarei, New Zealand (Lon:174.35 Lat:-35.76) to Gibraltar, Gibraltar (Lon: -5.35 Lat:36.13):

distVincentyEllipsoid(c(174.35,-35.76),c(-5.35,36.13))

> NA

...but:

Valetta, Malta (Lon:14.51 Lat:35.89) to Gibraltar, Gibraltar (Lon: -5.35 Lat:36.13):

distVincentyEllipsoid(c(14.51,35.89),c(-5.35,36.13))

> 1787499

...and:

Whangarei, New Zealand (Lon:174.35 Lat:-35.76) to Valetta, Malta (Lon:14.51 Lat:35.89):

distVincentyEllipsoid(c(174.35,-35.76),c(14.51,35.89))

> 18207301

(Both of these outputs are correct according to Google Maps - note these values are in metres.)

Any help would be greatly appreciated.

EDIT: Problem solved – Vincenty's distance method fails for near-antipodal points.

www
  • 38,575
  • 12
  • 48
  • 84
Anthony S.
  • 361
  • 1
  • 11
  • 1
    Puzzling behavior from what most people would call a seldom-used function might be better directed to the package maintainer: Try: `maintainer("geosphere")` – IRTFM Apr 06 '17 at 17:59
  • 1
    And especially weird since `distHaversine(c(174.35,-35.76),c(-5.35,36.13))` yields an answer. – SymbolixAU Apr 06 '17 at 21:52
  • Good point @42, I'll get in touch. – Anthony S. Apr 11 '17 at 14:18
  • And thanks @SymbolixAU for pointing that out. Using Haversine might be a good enough solution until my issue with the Vincenty function is resolved. – Anthony S. Apr 11 '17 at 14:18
  • 1
    After a bit more research, it appears that Vincenty's method can fail for (near) anitpodal points, and [this website](http://www.movable-type.co.uk/scripts/latlong-vincenty.html) suggests it can fail with distances greater than 19,936 km. Which makes sense in your example as the Haversine distance is 19,988 km – SymbolixAU Apr 12 '17 at 03:08
  • Oh, wow – it hadn't even occurred to me that it may be an issue with Vincenty's method itself. Thanks for looking into it! That seems to resolve my quandary. – Anthony S. Apr 12 '17 at 05:44
  • 1
    no worries; it was useful for me to know too for my [`spatial.data.table`](https://github.com/SymbolixAU/spatial.data.table) package I'm putting together (the [calculations](https://github.com/SymbolixAU/spatial.data.table/blob/master/R/SpatialCalculations.R) in particular). – SymbolixAU Apr 12 '17 at 05:51

1 Answers1

4

Use distGeo instead of distVincentyEllipsoid. This was added to geosphere in version 1.5 and it computes the geodesic distance using a better algorithm (by me!) than Vincenty's. This returns accurate results for all pairs of points. In particular, Whangarei, New Zealand (Lon:174.35 Lat:-35.76) to Gibraltar, Gibraltar (Lon: -5.35 Lat:36.13)

> distGeo(c(174.35,-35.76),c(-5.35,36.13))
[1] 19958627

This demonstration page shows examples of near antipodal geodesics displayed in Google Maps.

cffk
  • 1,839
  • 1
  • 12
  • 17