通过 R 中的 SpatialPolygonsDataFrame 对象列表提取栅格答案

【问题标题】：Extract raster by a list of SpatialPolygonsDataFrame objects in R通过 R 中的 SpatialPolygonsDataFrame 对象列表提取栅格
【发布时间】：2018-11-10 18:05:52
【问题描述】：

我正在尝试从单个大文件中为存储在列表中的 R 中的各种 SpatialPolygonsDataFrames (SPDF) 对象提取求和的栅格单元格值，然后将提取的值添加到 SPDF 对象属性表中。我想迭代这个过程，但不知道该怎么做。我找到了存储在单个 SPDF 对象中的多个多边形的有效解决方案（请参阅：https://gis.stackexchange.com/questions/130522/increasing-speed-of-crop-mask-extract-raster-by-many-polygons-in-r），但不知道如何将crop>mask>extract 过程应用于 SPDF 对象列表，每个包含多个多边形。这是一个可重现的示例：

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

以下是我想以 forloop 格式为整个列表运行的内容。对于列表中的单个 SPDF，以下一系列函数似乎有效：

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

但是当我尝试隔离 SPDF 列表 (raster::crop) 的第一个函数时，它不会返回光栅对象：

crop1 <- crop(c, extent(boundl[[1]])) #correctly returns object class 'RasterLayer'
cropl <- lapply(boundl, crop, c, extent(boundl)) #incorrectly returns objects of class 'SpatialPolygonsDataFrame'

当我尝试隔离 SPDF 列表 (raster::mask) 的掩码函数时，它返回错误：

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

我想纠正这些错误，并在一个循环中有效地迭代整个过程（即crop>mask>extract>将提取的值添加到 SPDF 属性表。我对 R 真的很陌生不知道从这里去哪里。请帮忙！

【问题讨论】：

bound1 未定义。

标签： r polygon spatial raster

【解决方案1】：

一种方法是采用有效的方法，只需将所需的“crop -> mask -> extract -> add”放入 for 循环中：

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
}

可以通过并行执行来加速循环，例如，使用 R 包 foreach。相反，使用lapply() 而不是 for 循环的速度增益会很小。

为什么会出现错误：

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

将函数crop() 应用于列表boundl 的每个元素。执行的操作是

tmp <- crop(boundl[[1]], c)
## test if equal to first element
all.equal(cropl[[1]], tmp) 
[1] TRUE

要获得所需的结果，请使用

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

注意：

使用c 表示一个R 对象是一个不好的选择，因为它很容易与c() 混淆。

【讨论】：