提取包含在一个SpatialPolygons中的栅格单元格数量

Question

我将创建一个函数，用于计算SpatialPolygonsDataframe对象的多边形内的栅格单元格数量，并将该值作为新列添加到其中，而不使用循环。 我找不到如何做到这一点...

以下是我的代码：

library(sp)
library(raster)
# Create a SpatialPolygonsDataframe and a raster objets to overlay
# Polygons
p1 <- rbind(c(-180,-20), c(-140,55), c(-50, 0), c(-140,-60), c(-180,-20))
p2 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0))
p3 <- rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0))
polys <- as(spPolygons(p1, p2, p3), "SpatialPolygonsDataFrame")
# Raster
p4 <- rbind(c(-180,10), c(0,90), c(40,90), c(145,-10),  
        c(-25, -15), c(-180,0), c(-180,10))
rpol <- spPolygons(p4)
r <- raster(ncol=90, nrow=45)
grd <- rasterize(rpol, r, fun=sum)

# Function to count the share of occupied cells per polygon
myFunction <- function(polys,grd){
    over <- crop(mask(grd, polys), polys)
    share <- length(over[!is.na(over)]) / ncell(over) * 100
  return(share)
}  

myFunction(polys,grd) 

非常感谢您的帮助！

- Marion LT

Answer 1

您可以使用函数raster::extract。默认情况下，它会提取SpatialPolygonDataFrame中每个多边形的单元格信息。输出结果是一个列表，因此您可以计算列表中对象的数量。

library(sp)
library(raster)
# Create a SpatialPolygonsDataframe and a raster objets to overlay
# Polygons
p1 <- rbind(c(-180,-20), c(-140,55), c(-50, 0), c(-140,-60), c(-180,-20))
p2 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0))
p3 <- rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0))
polys <- as(spPolygons(p1, p2, p3), "SpatialPolygonsDataFrame")

# Create raster
r <- raster(ncol=90, nrow=45)
values(r) <- 1:ncell(r)
# Extract values
r.ext <- extract(r, polys)
# Count number of cell for each polygon
lengths(r.ext)

- Sébastien Rochette

Answer 2

值得注意的是，使用terra可以大大加快这个过程。在我的情况下，运行类似这样的代码大约快80倍：

p1 <- rbind(c(-180,-20), c(-140,55), c(-50, 0), c(-140,-60), c(-180,-20))
p2 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0))
p3 <- rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0))
polys <- as(spPolygons(p1, p2, p3), "SpatialPolygonsDataFrame")
r <- raster(ncol=90, nrow=45)
values(r) <- 1:ncell(r)

v <- terra::vect(polys)
r <- terra::rast(r)
cells_per_poly <- terra::extract(r, v) %>% group_by(ID) %>% summarise(n=n())

- smartse