Search code examples
rdistanceraster

Does gridDistance calculate the geodetic distance when working in LatLong? (raster package in R)

gridDistance function in R (raster package) produce a map of distances from cell with a given value. The function gives the distance in the units of the current projection, or in meters if working in LatLong. I wonder if using LatLong the distance is "corrected" (geodetic). Otherwise I probably need to re-project to an equidistant projection. Thanks

Solution

  • Looking into the source of the truly excellent raster package it looks like your intuition is correct...

    # From gridDistance() when lonlat is true
if (lonlat) {
    distance <- pointDistance(

# From pointDistance()
if (! longlat ) {
    return( .planedist(p1[,1], p1[,2], p2[,1], p2[,2]) )
} else { 
    return( .haversine(p1[,1], p1[,2], p2[,1], p2[,2], r=6378137) )
}

# Finally from .haversine()
.haversine <- function(x1, y1, x2, y2, r=6378137) {
adj <- pi / 180
x1 <- x1 * adj
y1 <- y1 * adj
x2 <- x2 * adj
y2 <- y2 * adj
x <- sqrt((cos(y2) * sin(x1-x2))^2 + (cos(y1) * sin(y2) - sin(y1) * cos(y2) * cos(x1-x2))^2)
y <- sin(y1) * sin(y2) + cos(y1) * cos(y2) * cos(x1-x2)
return ( r * atan2(x, y) )
}

    So in short, yes, you are getting super accurate geodetic distances.