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

Question

我正在尝试从单个大文件中提取 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) 隔离第一个函数时，它不会 return 光栅对象:

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) 隔离掩码函数时，它 return 出现错误：

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 的新手，不知道从这里去哪里。请帮忙！

Answer 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() 混淆。

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

Extract raster by a list of SpatialPolygonsDataFrame objects in R

r

polygon

raster

spatial