# Missing 1.6% of expected area after buffering SpatialPoints from vector of widths

Goal: I am buffering a random set of points (i.e. centroids) with different widths from a vector of predefined areas to make a random null model of polygons in a given extent. I am iterating through a 1000 object list go generate the distribution, but have provided an example for only one object since my issue occurs at that level.

Issue: The total area of the buffered centroids is 32.19617 hectares less than the original total area.

Question: How is this area being cut? I feel like I am missing something about how `gBuffer()` works. I thought it might be a planar vs geodesic issue, but read that `regeos` will only handle projected spatial objects (which I have). It doesn't seem like a rounding issue because the difference is spread across 101 points (0.318774 hectares per point). Any suggestions for how else to troubleshoot?

The coordinate system used is NAD 1983 Albers. The PROJ.4 projection string is:
"+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0

Example

``````library(rgdal)
library(rgeos)

#the predefined areas
areas <- c(21.92466, 0.13743, 22.4576, 3.10855, 13.33003, 19.6971, 8.65893,
1.36051, 1.97611, 0.4128, 5.09534, 13.08232, 0.65519, 6.69976,
1.31828, 61.96348, 2.55861, 98.36752, 17.00504, 35.09418, 63.48043,
15.28773, 9.96236, 0.75299, 7.92188, 9.88541, 1.44117, 9.31345,
0.83406, 113.84276, 4.55541, 4.61159, 4.3415, 0.78915, 11.73815,
0.59766, 0.94186, 9.81746, 30.88698, 3.93619, 1.50583, 39.07577,
33.18197, 39.31604, 29.86065, 6.17087, 66.89905, 4.66883, 64.63708,
2.83463, 38.12346, 14.0179, 25.69435, 26.97518, 26.86608, 8.11912,
32.90833, 47.68297, 7.36192, 0.8557, 1.35078, 33.75896, 4.10684,
1.33971, 26.18847, 1.68767, 34.18756, 18.59203, 4.13611, 23.08783,
27.20771, 0.63805, 16.71614, 10.73375, 15.9914, 61.96926, 8.12024,
1.26379, 14.14531, 11.27807, 9.35613, 5.94766, 92.88387, 17.1014,
4.05735, 15.91729, 10.92878, 3.23904, 4.15417, 59.53967, 0.41474,
6.3664, 16.00491, 0.72044, 14.36096, 1.10064, 2.52104, 6.35359,
33.68302, 176.33294, 12.89539)

extent <- readOGR(dsn = "SanDiegoCounty.shp", layer = "SanDiegoCounty")

#sample points within the extent equal to the total number in 'areas'
centroids <- spsample(x = extent, n = length(areas), type = "random")

#buffer centroids with width derived from predefined areas vector
buffered.centroids <- gBuffer(spgeom = centroids, byid = TRUE, id = names(centroids), width = sqrt(areas/pi))

#the total predefined area: 1966.97647
sum(areas)

#the total area of buffered centroids: 1934.78029709348
buffered.centroid.area <- sum(sapply(X = buffered.centroids@polygons, FUN = function(z){sum(z@area)}))

#difference in total areas: 32.1961729065185
sum(areas)-buffered.centroid.area
``````

Link to `extent` shapefile on Google Drive: extent.zip

• Please Edit the question to specify the units of the areas, and the coordinate system (.prj) of the data file. Nov 17, 2017 at 21:52
• Possibly different vertices that make up the boundaries. They can affect the area calculations. Nov 17, 2017 at 23:32

The answer by @obrl_soil is spot on. If you want to increase `gBuffer`'s approximation of circles, you can increase `quadsegs`. Here is a table of these errors:

``````library(rgeos)
pe <- c()
for(qs in 1:20) {
pe[qs] <- 100 * (1.0 - gArea(gBuffer(p0, quadsegs=qs)) / pi)
}

# Shows these messages ... and plot
``````

Two things you can see here is that the error for quadsegs=5 is 1.64%, and that the error levels off near 0.1% at quadsegs=20. You can push this to a crazy quadsegs=1000 to get 0.0000411%, if you really want, but you can't get exactly 0% error that you would expect from perfect circle buffers.

• very useful addition! I've added some comparative speed tests as well. Nov 21, 2017 at 3:40

gBuffer on points won't generate a true circle by default, but rather a regular 20-sided polygon (the setting 'quadsegs' is relevant here and defaults to 5). As the wiki article for Icosagon points out, the area of such a shape is very slightly smaller than the true circle that contains it.

If its important to get closer to a perfect circle, try increasing quadsegs to 10 or even 20, but bear in mind that performance will be affected:

``````library(sp)
library(rgeos)
library(sf)
library(microbenchmark)
library(ggplot2)

# 1000 random points
pts <- matrix(c(replicate(2, sample(10000, size = 1000, replace = FALSE))),
ncol = 2, byrow = F)
pts <- st_multipoint(pts)
pts <- st_sfc(pts)
pts <- st_cast(pts, 'POINT')
pts_sp <- as(pts, 'Spatial')

# 1. test speed as quadsegs increases in gBuffer
test_rg <- microbenchmark(
"5"  = gBuffer(pts_sp),
"10" = gBuffer(pts_sp, quadsegs = 10),
"20" = gBuffer(pts_sp, quadsegs = 20),
"50" = gBuffer(pts_sp, quadsegs = 50)
)

# 2. test equivalent function in sf package
test_sf <- microbenchmark(
"5"  = st_buffer(pts, dist = 1, nQuadSegs = 5),
"10" = st_buffer(pts, dist = 1, nQuadSegs = 10),
"20" = st_buffer(pts, dist = 1, nQuadSegs = 20),
"50" = st_buffer(pts, dist = 1, nQuadSegs = 50)
)

ggplot() +
stat_summary(data = test_rg,
aes(x = expr, y = time/1000000, group = 1, colour = 'red'),
size = 1.2, fun.y = 'median', geom = 'line') +
stat_summary(data = test_sf,
aes(x = expr, y = time/1000000, group = 1, colour = 'blue'),
size = 1.2, fun.y = 'median', geom = 'line') +
labs(x = 'n quadsegs', y = 'milliseconds') +
ggtitle('Buffer function speed comparison',
subtitle = '1000 points, 100 runs, buffer dist 1') +
scale_colour_manual(name = 'package::function',
values = c('blue'='blue','red'='red'),
labels = c('sf::st_buffer','rgeos::gBuffer'),
guide = 'legend')
``````

sf and rgeos wrap the same C++ library so results will be identical. Note that rgeos::gBuffer defaults to 5 quadsegs, while sf::st_buffer defaults to 30, so the latter package sacrifices some speed for accuracy by default. sf's inherent speed advantage offsets that somewhat.