CODEWITHSUNDEEP
rdistancelatitude-longitude
Christoph
Christoph

Reputation: 7063

Why does geoshere give the wrong distance?

Just checkout the following example:

lat <- c(47.5, 47.5, 47.501)
lng <- c(9.3, 9.301, 9.3)
geosphere::distVincentyEllipsoid(lng[c(1,2)], lat[c(1,2)])
# [1] 5551424 Should be about 5000 m?
geosphere::distVincentyEllipsoid(lng[c(1,3)], lat[c(1,3)])
#[1] 5551583 Should be about 5000 m?

m <- leaflet() %>%
  addTiles(group = "OSM") %>%
  addProviderTiles("Stamen.TonerLite") %>%
  addLayersControl(
    baseGroups = c("OSM", "Stamen.TonerLite")) %>% 
  addCircleMarkers(lat = lat,
                   lng = lng,
                   color = "blue",
                   stroke = FALSE,
                   radius = 3,
                   fillOpacity = 0.7)
print(m)

To me it seem, the result is quite off although ????distVincentyEllipsoid state:

Distance value in the same units as the ellipsoid (default is meters)

Edit: Just as a simple crosscheck from here:

R = 6371  # radius of the earth in km
x = (lng[2] - lng[1]) * cos( 0.5*(lat[2]+lat[1]) )
y = lat[2] - lat[1]
d = R * sqrt( x*x + y*y ) # in km

Upvotes: 2

Views: 179

Answers (1)

Jaap
Jaap

Reputation: 83215

You are supplying distVincentyEllipsoid the wrong information. It should have the following format: distVincentyEllipsoid(p1, p2). From the help-file:

p1  longitude/latitude of point(s), in degrees 1; can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints* object
p2  as above

What you are doing is providing the two lng values to p1 and the two lat values to p2, whereas both points should be a lng/lat combination.

If you do:

distVincentyEllipsoid(c(lng[1], lat[1]), c(lng[2], lat[2]))

you will get:

[1] 75.34357

Upvotes: 3

Related Questions