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.
R
can do that--it's a great suggestion. But since for a spherical Earth model this is such a simple projection, it's simple enough to write the code directly.