2

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)

#don't know how to code SpatialPolygonsDataFrame for reproduciblity, but it can be downloaded at link below
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. – Vince Nov 17 '17 at 21:52
  • Possibly different vertices that make up the boundaries. They can affect the area calculations. – mkennedy Nov 17 '17 at 23:32
2

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)
p0 <- readWKT("POINT(0 0)")
pe <- c()
for(qs in 1:20) {
    pe[qs] <- 100 * (1.0 - gArea(gBuffer(p0, quadsegs=qs)) / pi)
    message(sprintf("quadsegs=%d : error=%.3g%%", qs, pe[qs]))
}
plot(1:20, pe, "b", xlab="quadsegs")

# Shows these messages ... and plot
quadsegs=1 : error=36.3%
quadsegs=2 : error=9.97%
quadsegs=3 : error=4.51%
quadsegs=4 : error=2.55%
quadsegs=5 : error=1.64%
quadsegs=6 : error=1.14%
quadsegs=7 : error=0.837%
quadsegs=8 : error=0.641%
quadsegs=9 : error=0.507%
quadsegs=10 : error=0.411%
quadsegs=11 : error=0.34%
quadsegs=12 : error=0.285%
quadsegs=13 : error=0.243%
quadsegs=14 : error=0.21%
quadsegs=15 : error=0.183%
quadsegs=16 : error=0.161%
quadsegs=17 : error=0.142%
quadsegs=18 : error=0.127%
quadsegs=19 : error=0.114%
quadsegs=20 : error=0.103%

enter image description here

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. – obrl_soil Nov 21 '17 at 3:40
2

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')

enter image description here

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.