Extract raster information from a spdf and for different buffer size in a loop

Question

For a project in geoprocessing, I need to extract information from a raster but for different buffer sizes. I have a raster of the population (100 m size pixel) and a SpatialPolygonDataFrame with the localisation of volcanoes. According to the volcanic explosivity index (VEI) of those volcanoes, I created a new column in the spdf with different buffer sizes (VEI = 1 --> buffer_size = 5000; VEI = 2 --> buffer_size = 10000, ...). Now, I want to extract the sum of pixel included inside the buffer. I want to do it the more efficient ! Then, I wanted to make a loop, but it doesn't work... Can somebody tell me why?

When I want to extract for a single buffer value, there is no problem (see the following code) :

people_number <- extract(x = population_proj, y = locations_spdf_proj, buffer = 10000, fun = sum)

But When I want to do that in a loop, I have a warning message (Error in .cellValues(x, y, ...)) :

  for (i in 1:length(locations_spdf_proj)){

  people[[i]] <- extract(x = population_proj,                            # raster layer
                         y = locations_spdf_proj[[i]],                   # SPDF with locations for buffer
                         buffer = locations_spdf_proj$buffer_size[[i]],  # buffer size
                         fun = sum)                                       # value to extract (sum of all pixels)

}

Answer 1

The buffer argument in raster::extract is intended for point feature classes, not polygons. You are going to need to create a new polygon feature class, using rgeod::gBuffer, and extract for the new set of buffered polygons.

The buffered results should remain ordered to the original polygons so, can just be merged later. If need to do this over multiple columns of "buffer by" values, you can just pipe the extracted sums into an empty list object and use do.call with cbind to merge the results into a data.frame that could be added back to the original polygons. For one set of variable buffers, you do not even need a loop to do this.

Please note that in the example I am using your object names for clarity sake but, please, in the future, shorten your naming convention as it is quite cumbersome and unnecessary.

• population_proj - raster data
• locations_spdf_proj - polygon data
• locations_spdf_proj$buffer_size - buffer size attribute
• locations_spdf_proj.buff - Buffered polygon data based on buffer_size
• locations_spdf_proj$pop.sum - New raster sum attribute in original polygons

First, add packages and create some example data.

library(raster)
library(sp)

locations_spdf_proj <- raster(res=1000, xmn=571823.6, xmx=616763.6, ymn=4423540, 
                     ymx=4453690, crs = CRS("+proj=utm +zone=12 +datum=NAD83 
                     +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
  population_proj <- locations_spdf_proj
  population_proj[] <- rpois(ncell(population_proj), lambda=1)  
  locations_spdf_proj[] <- runif(ncell(locations_spdf_proj)) * 10
  locations_spdf_proj <- rasterToPolygons(locations_spdf_proj, fun=function(x){x > 9.9})

Add some variable buffer distances to the polygons, buffer them and plot the results of the buffered and original polygons on top of the raster.

locations_spdf_proj$buffer_size <- runif(nrow(locations_spdf_proj), 1000,2500)      
locations_spdf_proj.buff <- rgeos::gBuffer(locations_spdf_proj, byid=TRUE,
                                    width=locations_spdf_proj$buffer_size)


plot(population_proj)
  plot(locations_spdf_proj.buff, add=TRUE)
  plot(locations_spdf_proj, add=TRUE)

Now, we can use extract to get the sums the buffered polygon and add a new column to the original polygon object.

( locations_spdf_proj$pop.sum <- extract( population_proj,
                      locations_spdf_proj.buff, fun=sum) )

  spplot(locations_spdf_proj, "pop.sum")