Nrthugu

I am trying to extract summed raster cell values from a single big file for various SpatialPolygonsDataFrames (SPDF) objects in R stored in a list, then add the extracted values to the SPDF objects attribute tables. I would like to iterate this process, and have no idea how to do so. I have found an efficient solution for multiple polygons stored in a single SPDF object (see: https://gis.stackexchange.com/questions/130522/increasing-speed-of-crop-mask-extract-raster-by-many-polygons-in-r), but do not know how to apply the crop>mask>extract procedure to a LIST of SPDF objects, each containing multiple polygons. Here is a reproducible example:

library(maptools)  ## For wrld_simpl

library(raster)



## Example SpatialPolygonsDataFrame

data(wrld_simpl) #polygon of world countries

bound <- wrld_simpl[1:25,] #country subset 1  

bound2 <- wrld_simpl[26:36,] #subset 2



## Example RasterLayer

c <- raster(nrow=2e3, ncol=2e3, crs=proj4string(wrld_simpl), xmn=-180, 

xmx=180, ymn=-90, ymx=90)

c <- 1:length(c)



#plot, so you can see it

plot(c)    

plot(bound, add=TRUE) 

plot(bound2, add=TRUE, col=3) 



#make list of two SPDF objects

boundl<-list()

boundl[[1]]<-bound1

boundl[[2]]<-bound2



#confirm creation of SPDF list

boundl

The following is what I would like to run for the entire list, in a forloop format. For a single SPDF from the list, the following series of functions seem to work:

clip1 <- crop(c, extent(boundl[[1]])) #crops the raster to the extent of the polygon, I do this first because it speeds the mask up

clip2 <- mask(clip1, boundl[[1]]) #crops the raster to the polygon boundary 

extract_clip <- extract(clip2, boundl[[1]], fun=sum)  

#add column + extracted raster values to polygon dataframe

boundl[[1]]@data["newcolumn"] = extract_clip

But when I try to isolate the first function for the SPDF list (raster::crop), it does not return a raster object:

crop1 <- crop(c, extent(boundl[[1]])) #correctly returns object class 'RasterLayer'

cropl <- lapply(boundl, crop, c, extent(boundl)) #incorrectly returns objects of class 'SpatialPolygonsDataFrame'

When I try to isolate the mask function for the SPDF list (raster::mask), it returns an error:

maskl <- lapply(boundl, mask, c) 

#Error in (function (classes, fdef, mtable)  : unable to find an inherited method for function ‘mask’ for signature ‘"SpatialPolygonsDataFrame", "RasterLayer"’

I would like to correct these errors, and efficiently iterate the entire procedure within a single loop (i.e., crop>mask>extract>add extracted values to SPDF attribute tables. I am really new to R and don't know where to go from here. Please help!

One approach is to take what is working and simply put the desired "crop -> mask -> extract -> add" into a for loop:

for(i in seq_along(boundl)) {

    clip1 <- crop(c, extent(boundl[[i]])) 

    clip2 <- mask(clip1, boundl[[i]])

    extract_clip <- extract(clip2, boundl[[i]], fun=sum)  

    boundl[[i]]@data["newcolumn"] <- extract_clip

}

One can speed-up the loop with parallel execution, e.g., with the R package foreach. Conversely, the speed gain of using lapply() instead of the for loop will be small.

Why the error occurs:

cropl <- lapply(boundl, crop, c, extent(boundl)) 

applies the function crop() to each element of the list boundl. The performed operation is

tmp <- crop(boundl[[1]], c)

## test if equal to first element

all.equal(cropl[[1]], tmp) 

[1] TRUE

To get the desired result use

cropl <- lapply(boundl, function(x, c) crop(c, extent(x)), c=c)

## test if the first element is as expected

all.equal(cropl[[1]], crop(c, extent(boundl[[1]]))) 

[1] TRUE

Note:

Using c to denote an R object is a bade choice, because it can be easily confused with c().