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