Euclidean and Geodesic Buffering in R

Question

In Understanding Geodesic Buffering, The Esri Geoprocessing Development Team distinguish between Euclidean and Geodesic Buffering. They conclude with "Euclidean buffering performed on projected feature classes can produce misleading and technically incorrect buffers. However, geodesic buffering will always produce results that are geographically accurate because geodesic buffers are not affected by the distortions introduced by projected coordinate systems".

I have to work with a point global dataset and the coordinates are unprojected (+proj=longlat +ellps=WGS84 +datum=WGS84). Is there a function to create a geodesic buffer in R when width is given in metric units? I am aware of gBuffer from rgeos package. This function creates a buffer in units of the spatial object that is used (example), so, I have to project the coordinates to be able to create a buffer of desired X km. Projecting and then applying a gBuffer means actually making an Euclidean buffer as opposed to a Geodesic one that I need. Below is some code to illustrate my concerns:

require(rgeos)
require(sp)
require(plotKML)

# Generate a random grid-points for a (almost) global bounding box
b.box <- as(raster::extent(120, -120, -60, 60), "SpatialPolygons")
proj4string(b.box) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
set.seed(2017)
pts <- sp::spsample(b.box, n=100, type="regular")
plot(pts@coords)

# Project to Mollweide to be able to apply buffer with `gBuffer` 
# (one could use other projection)
pts.moll <- sp::spTransform(pts, CRSobj = "+proj=moll")
# create 1000 km buffers around the points
buf1000km.moll <- rgeos::gBuffer(spgeom = pts.moll, byid = TRUE, width = 10^6)
plot(buf1000km.moll)
# convert back to WGS84 unprojected
buf1000km.WGS84 <- sp::spTransform(buf1000km.moll, CRSobj = proj4string(pts))
plot(buf1000km.WGS84) # distorsions are present
# save as KML to better visualize distorted Euclidian buffers on Google Earth
plotKML::kml(buf1000km.WGS84, file.name = "buf1000km.WGS84.kml")

The image below depicts the distorted Euclidian buffers (1000 km radius) produced with the code from above.

Robert J. Hijmans in Introduction to the ”geosphere” package, section 4 Point at distance and bearing gives an example of how to make "circular polygons with a fixed radius, but in longitude/latitude coordinates", which I think can be called a "geodesic buffer". I fallowed this idea and I wrote some code that hopefully does the right thing, but I wonder if there is already some geodesic-buffer R function in some package that allows metric radius as input:

require(geosphere)

make_GeodesicBuffer <- function(pts, width) {
    ### A) Construct buffers as points at given distance and bearing
    # a vector of bearings (fallows a circle)
    dg <- seq(from = 0, to = 360, by = 5)

    # Construct equidistant points defining circle shapes (the "buffer points")
    buff.XY <- geosphere::destPoint(p = pts, 
                                    b = rep(dg, each = length(pts)), 
                                    d = width)

    ### B) Make SpatialPolygons
    # group (split) "buffer points" by id
    buff.XY <- as.data.frame(buff.XY)
    id  <- rep(1:length(pts), times = length(dg))
    lst <- split(buff.XY, id)

    # Make SpatialPolygons out of the list of coordinates
    poly   <- lapply(lst, sp::Polygon, hole = FALSE)
    polys  <- lapply(list(poly), sp::Polygons, ID = NA)
    spolys <- sp::SpatialPolygons(Srl = polys, 
                                  proj4string = CRS(as.character("+proj=longlat +ellps=WGS84 +datum=WGS84")))
    # Disaggregate (split in unique polygons)
    spolys <- sp::disaggregate(spolys)
    return(spolys)
}

buf1000km.geodesic <- make_GeodesicBuffer(pts, width=10^6)
# save as KML to visualize geodesic buffers on Google Earth
plotKML::kml(buf1000km.geodesic, file.name = "buf1000km.geodesic.kml")

The image below depicts the Geodesic buffers (1000 km radius).

Edit 2019-02-12: For convenience, I wrapped a version of the function in the geobuffer package. Feel free to contribute with pull requests.

Answer 1

For most purposes it will be accurate enough to use a spherical model of the earth--and the coding will be easier and the calculations much faster.

Following suggestions by M. Kennedy in a comment, this solution buffers the North Pole (which is easy: the buffer boundary lies at a fixed latitude) and then rotates this buffer into any desired location.

The rotation is effected by converting the original buffer to geocentric Cartesian (XYZ) coordinates, rotating those with a (fast) matrix multiplication along the Prime Meridian to the target latitude, converting its coordinates back to Geographic (lat-lon), and then spinning the buffer around the Earth's axis simply by adding the target longitude to each second coordinate.

Why do it in two steps when (ordinarily) a single matrix multiplication would work? Because there is no need for special code to identify the breaks at the +/-180 degree meridian. Instead, this approach can generate longitudes beyond the original range (whether -180 to 180 degrees or 0 to 360 or whatever), but by so doing, standard polygon-drawing procedures (and other analytical procedures) will work fine without modification. If you must have longitudes in a given range, simply reduce them modulo 360 degrees at the very end: that's fast and easy.

When buffering points, typically all buffers have the same radius. This modular solution allows for a speedup in this case: you can buffer the North Pole and then convert it into XYZ coordinates once and for all. Buffering each point thereby requires on (very fast) matrix multiplication, conversion back to lat-lon coordinates, and shifting of the longitudes (also very fast). Expect to generate about 10,000 high-resolution buffers (360 vertices) per second.

This R code shows the details. Since its purpose is illustration, it uses no add-on packages; nothing is hidden or buried. It includes a test in which a set of random points is generated, buffered, and displayed using its raw lat-lon (Geographic) coordinates. Here's the output:

degrees.to.radians <- function(phi) phi * (pi / 180)
radians.to.degrees <- function(phi) phi * (180 / pi)
#
# Create a 3X3 matrix to rotate the North Pole to latitude `phi`, longitude 0.
# Solution: A rotation is a linear map, and therefore is determined by its
#           effect on a basis.  This rotation does the following:
#           (0,0,1) -> (cos(phi), 0, sin(phi))  {North Pole (Z-axis)}
#           (0,1,0) -> (0,1,0)                  {Y-axis is fixed}
#           (1,0,0) -> (sin(phi), 0, -cos(phi)) {Destination of X-axis}
#
rotation.create <- function(phi, is.radians=FALSE) {
  if (!is.radians) phi <- degrees.to.radians(phi)
  cos.phi <- cos(phi)
  sin.phi <- sin(phi)
  matrix(c(sin.phi, 0, -cos.phi, 0, 1, 0, cos.phi, 0, sin.phi), 3)
}
#
# Convert between geocentric and spherical coordinates for a unit sphere.
# Assumes `latlon` in degrees.  It may be a 2-vector or a 2-row matrix.
# Returns an array with three rows for x,y,z.
#
latlon.to.xyz <- function(latlon) {
  latlon <- degrees.to.radians(latlon)
  latlon <- matrix(latlon, nrow=2)
  cos.phi <- cos(latlon[1,])
  sin.phi <- sin(latlon[1,])
  cos.lambda <- cos(latlon[2,])
  sin.lambda <- sin(latlon[2,])
  rbind(x = cos.phi * cos.lambda,
        y = cos.phi * sin.lambda,
        z = sin.phi)
}
xyz.to.latlon <- function(xyz) {
  xyz <- matrix(xyz, nrow=3) 
  radians.to.degrees(rbind(phi=atan2(xyz[3,], sqrt(xyz[1,]^2 + xyz[2,]^2)),
                           lambda=atan2(xyz[2,], xyz[1,])))
}
#
# Create a circle of radius `r` centered at the North Pole, oriented positively.
# `r` is measured relative to the sphere's radius `R`.  For the authalic Earth,
# r==1 corresponds to 6,371,007.2 meters.
#
# `resolution` is the number of vertices to use in a polygonal approximation.
# The first and last vertex will coincide.
#
circle.create <- function(r, resolution=360, R=6371007.2) {
  phi <- pi/2 - r / R                       # Constant latitude of the circle
  resolution <- max(1, ceiling(resolution)) # Assures a positive integer
  radians.to.degrees(rbind(phi=rep(phi, resolution+1),
                           lambda=seq(0, 2*pi, length.out = resolution+1)))
}
#
# Rotate around the y-axis, sending the North Pole to `phi`; then
# rotate around the new North Pole by `lambda`.
# Output is in geographic (spherical) coordinates, but input points may be
# in Earth-centered Cartesian or geographic.
# No effort is made to clamp longitudes to a 360 degree range.  This can 
# facilitate later computations.  Clamping is easily done afterwards if needed:
# reduce the longitude modulo 360 degrees.
#
rotate <- function(p, phi, lambda, is.geographic=FALSE) {
  if (is.geographic) p <- latlon.to.xyz(p)
  a <- rotation.create(phi)   # First rotation matrix
  q <- xyz.to.latlon(a %*% p) # Rotate the XYZ coordinates
  q + c(0, lambda)            # Second rotation
}
#------------------------------------------------------------------------------#
#
# Test.
#
n <- 50                  # Number of circles
radius <- 1.5e6          # Radii, in meters
resolution <- 360
set.seed(17)             # Makes this code reproducible

#-- Generate random points for testing.
centers <- rbind(phi=(rbeta(n, 2, 2) - 1/2) * 180,
                 lambda=runif(n, 0, 360))

system.time({
  #-- Buffer the North Pole and convert to XYZ once and for all.
  p.0 <- circle.create(radius, resolution=resolution) 
  p <- latlon.to.xyz(p.0)

  #-- Buffer the set of points (`centers`).
  circles <- apply(centers, 2, function(center) 
    rotate(p, center[1], center[2]))

  #-- Convert into an array indexed by (latlon, vertex, point id).
  circles <- array(circles, c(2, resolution+1, n))
})
#
# Display the buffers (if there are not too many).
#
if (n <= 1000) {
  #-- Create a background map area and graticule.
  xlim <- range(circles[2,,]) # Extent of all longitudes in the buffers
  plot(xlim, c(-90, 90), type="n", xlim=xlim, ylim=c(-90,90), asp=1,
       xlab="Longitude", ylab="Latitude",
       main=paste(n, "Random Disks of Radius", signif(radius/1e3, 3), "Km"),
       sub="Centers shown with gray dots")
  abline(v=seq(-360, 720, by=45), lty=1, col="#d0d0d0")
  abline(h=seq(-90, 90, by=30), lty=1, col="#d0d0d0")

  #-- Display the buffers themselves.
  colors <- terrain.colors(n, alpha=1/3) # Vary their colors
  invisible(sapply(1:n, function(i) {
    polygon(circles[2,,i], circles[1,,i], col=colors[i])
  }))

  #-- Show the original points (and, optionally, labels).
  points(centers[2,], centers[1,], pch=21, bg="Gray", cex=min(1, sqrt(25/n)))
  # text(centers[2,], centers[1,], labels=1:n, cex=min(1, sqrt(100/n)))
}