由一个或多个多边形覆盖的栅格单元的一部分:是否有更快的方法(在R中)？

Question

图片比文字好,所以请看看

我拥有的是什么

• 一个RasterLayer对象(此处仅填充随机值仅用于说明目的,实际值无关紧要)
• 一个SpatialPolygons对象,里面有很多很多的多边形

您可以使用以下代码重新创建我用于图像的示例数据:

library(sp)
library(raster)
library(rgeos)

# create example raster
r <- raster(nrows=10, ncol=15, xmn=0, ymn=0)
values(r) <- sample(x=1:1000, size=150)

# create example (Spatial) Polygons
p1 <- Polygon(coords=matrix(c(50, 100, 100, 50, 50, 15, 15, 35, 35, 15), nrow=5, ncol=2), hole=FALSE)
p2 <- Polygon(coords=matrix(c(77, 123, 111, 77, 43, 57, 66, 43), nrow=4, ncol=2), hole=FALSE)
p3 <- Polygon(coords=matrix(c(110, 125, 125, 110, 67, 75, 80, 67), nrow=4, ncol=2), hole=FALSE)

lots.of.polygons <- SpatialPolygons(list(Polygons(list(p1, p2, p3), 1)))
crs(lots.of.polygons) <- crs(r) # copy crs from raster to polygons (please ignore any potential problems related to projections etc. for now)


# plot both
plot(r) #values in this raster for illustration purposes only
plot(lots.of.polygons, add=TRUE)

对于栅格中的每个单元格,我想知道一个或多个多边形覆盖了多少单元格.或者实际上:栅格单元格内所有多边形的区域,没有相关单元格之外的区域.如果有多个多边形重叠一个单元格,我只需要它们的组合区域.

以下代码执行我想要的操作,但使用实际数据集运行需要一周多的时间:

# empty the example raster (don't need the values):
values(r) <- NA

# copy of r that will hold the results
r.results <- r

for (i in 1:ncell(r)){
  r.cell <- r # fresh copy of the empty raster
  r.cell[i] <- 1 # set the ith cell to 1
  p <- rasterToPolygons(r.cell) # create a polygon that represents the i-th raster cell
  cropped.polygons <- gIntersection(p, lots.of.polygons) # intersection of i-th raster cell and all SpatialPolygons

  if (is.null(cropped.polygons)) {
    r.results[i] <- NA # if there's no polygon intersecting this raster cell, just return NA ...
  } else{
    r.results[i] <- gArea(cropped.polygons) # ... otherwise return the area
  }
}

plot(r.results)
plot(lots.of.polygons, add=TRUE)

我可以通过使用sapply而不是for-loop 来挤出更多的速度,但瓶颈似乎在其他地方.整个方法感觉很尴尬,我想知道我是否错过了一些明显的东西.起初我认为rasterize()应该能够轻松地做到这一点,但我无法弄清楚要在fun=论证中加入什么.有任何想法吗？

Answer 1

[编辑]

也许gIntersection(..., byid = T)with gUnaryUnion(lots.of.polygons)（它们使您能够一次处理所有细胞）比 for 循环更快（如果gUnaryUnion()花费太多时间，这是一个坏主意）。

r <- raster(nrows=10, ncol=15, xmn=0, ymn=0)
set.seed(1); values(r) <- sample(x=1:1000, size=150)
rr <- rasterToPolygons(r)

# joining intersecting polys and put all polys into single SpatialPolygons
lots.of.polygons <- gUnaryUnion(lots.of.polygons)   # in this example, it is unnecessary

gi <- gIntersection(rr, lots.of.polygons, byid = T)

ind <- as.numeric(do.call(rbind, strsplit(names(gi), " "))[,1])   # getting intersected rr's id
r[] <- NA
r[ind] <- sapply(gi@polygons, function(x) slot(x, 'area'))  # a bit faster than gArea(gi, byid = T)

plot(r)
plot(lots.of.polygons, add=TRUE)