Skip to content

Instantly share code, notes, and snippets.

@walkerke
Forked from valentinitnelav/make_GeodesicBuffer.R
Last active December 13, 2019 15:51
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save walkerke/4a7ec533f5de2f78d0305cf1956f73d2 to your computer and use it in GitHub Desktop.
Save walkerke/4a7ec533f5de2f78d0305cf1956f73d2 to your computer and use it in GitHub Desktop.
Geodesic buffer for point data
make_GeodesicBuffer <- function(XY.dg, dg.step=5, dst.m, crs){
#######################################################################################################
## Function to make a circle-buffer around given points (long-lat coordinates)
## Is different from rgeos::gBuffer() by the fact that it allows the user to create
## a geodesic buffer with a width expressed in metric units.
## Otherwise the user needs to project and apply an Euclidian buffer with rgeos::gBuffer(),
## which will introduce distortions that vary greatly with latitude and the radius of the circle buffer.
## See also my question addressed here:
## https://gis.stackexchange.com/questions/250389/euclidean-and-geodesic-buffering-in-r
##
## ARGUMENTS
## - XY.dg:
## matrix or spatial point object (SpatialPointsDataFrame) with long-lat coordinates
## (column 1 = longitude, column 2 = latitude)
## - dg.step:
## bearing (direction) step in degrees; dictates the point density of the circle-buffer edge
## - dst.m:
## distance (radious of circle buffer) in meters
## - crs:
## string of class CRS-class (CRS-class {sp})
## e.g. CRS(as.character("+proj=longlat +ellps=WGS84 +datum=WGS84"))
##
## RETURNS SpatialPolygons circle buffer
#######################################################################################################
# ---------------------------------------
# Check for valid input and packages
# ---------------------------------------
# Check for packages
if (!requireNamespace("geosphere", quietly = TRUE))
stop("Please install package {geosphere}")
if (!requireNamespace("data.table", quietly = TRUE))
stop("Please install package {data.table}")
if (!requireNamespace("sp", quietly = TRUE))
stop("Please install package {sp}")
# Check if XY.dg argument is of expected class
is.class.ok <- inherits(XY.dg, c("SpatialPoints","SpatialPointsDataFrame","matrix"))
if (!is.class.ok) stop("For XY.dg argument expecting 'SpatialPoints', 'SpatialPointsDataFrame' or 'matrix' class")
# depending on class of XY.dg, get number of points
if( inherits(XY.dg, c("SpatialPoints","SpatialPointsDataFrame")) ){
N.points <- length(XY.dg) # if is an sp object, take length
# Also if spatial object is projected then stop with error
if (is.projected(XY.dg)) stop("Spatial object XY.dg is projected; expects unprojected coordinates (long-lat)")
} else {
# if matrix - some routine checking of coordinates values
if (dim(XY.dg)[2] != 2) stop("Expecting XY.dg to be a 2 columns matrix (column1=longitude, column2=latitude)")
if ( ! (range(XY.dg[,1]) %between% c(-180, 180) && range(XY.dg[,2]) %between% c(-90, 90)) )
stop("Expects unprojected coordinates in XY.dg (longitude between -180 & 180, latitude between -90 & 90)")
N.points <- dim(XY.dg)[1] # if is a matrix, take number of rows
}
# ===========================================================
# A) Construct buffers as points at given distance and bearing
# ===========================================================
# a vector of bearings (fallows a circle)
dg <- seq(from = 0, to = 360, by = dg.step)
# # Construct equidistant points defining circle shapes (the "buffer points")
# Inspired from section "Point at distance and bearing" from
# Robert J. Hijmans in Introduction to the ”geosphere” package at:
# https://cran.r-project.org/web/packages/geosphere/vignettes/geosphere.pdf
buff.XY <- geosphere::destPoint(p = XY.dg,
b = rep(dg, each = N.points),
d = dst.m)
# ===========================================================
# B) Make SpatialPolygon from the points above
# ===========================================================
buff.XY <- data.table(buff.XY)
# group (split) "buffer points" by id
# add column which indicates to which point ID from N.points each circle-buffer point belongs to
buff.XY[, id := rep(1:N.points, times = length(dg))]
# buff.XY[, bearing := rep(dg, each = dim(XY.dg)[1])] # for debugging reasons
# split in a list of data tables with two columns: lon and lat
lst <- split(buff.XY[,.(lon,lat,id)], by = "id", keep.by = FALSE)
# 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(crs))
# Disaggregate (split in unique polygons)
spolys.buff <- sp::disaggregate(spolys)
return(spolys.buff)
}
@valentinitnelav
Copy link

Hi @walkerke,

Thanks for linking my function. Is currently wrapped in the geobuffer package.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment