使用一个shapefile（具有许多多边形）批处理/提取一个栅格的原始数据？

Question

如果我们有一个栅格，比如说一个国家的整数高程数据， 和一个多边形 shapefile，比如将那个国家的 300 个流域离散化，每个流域都有一个唯一的名称，我们如何最容易地为它们都获得这样的输出？

basinID, gridcellelev
a, 320
a, 321
a, 320
b, 17
b, 18
b, 19

最繁重的方法似乎是将单个 shapefile 过滤/转换为 300 个 shapefile， 将栅格 300 次裁剪为 300 个 uniqueID 栅格，将它们读回，为每个盆地生成单独的表格，然后将它们rbind()ing 在一起。

另一方面，理想的方法似乎是跳过文件生成，不保存 xy 数据，并仅使用一个栅格和一个 shapefile 创建同一个表 - 通过以某种方式选择盆地中的单元格，用盆地ID，创建一个表，丢失坐标，并不断迭代和附加该表，直到第 300 个盆地。

我不是在寻找任何统计数据，只是在某些标准剪辑方法中列出的网格单元标高的原始数据列表。我相信来自 ArcMap 的标准栅格剪辑属性表输出不是原始数据，而是像元的计数/频率。那也行。

我不知道最小限度地复制栅格和多边形 shapefile， 所以我很感激任何提示/库/函数/示例。作为起点：

library(tidyverse)
library(raster)
library(rgdal)
library(sf)

elev_raster <- raster("spain_elev_meters.tif") #integer raster
basins <- readOGR("spainbasins.shp", "spainbasins") %>% st_as_sf() #unique basin ID column: `basinID` 

除非建议使用其他起始格式。

非常感谢任何提示！

Answer 1

如果你搜索“提取光栅形状文件”，我想你会找到答案

示例数据：

library(raster)
p <- shapefile(system.file("external/lux.shp", package="raster"))
r <- raster(p)
values(r) <- 1:ncell(r)

p 有 12 个多边形

p
#class       : SpatialPolygonsDataFrame 
#features    : 12 
#extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
#crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
#variables   : 5
#names       : ID_1,     NAME_1, ID_2,   NAME_2, AREA 
#min values  :    1,   Diekirch,    1, Capellen,   76 
#max values  :    3, Luxembourg,   12,    Wiltz,  312

解决方案 1. 提取值并应用函数。例如mean

x <- extract(r, p, fun=mean, na.rm=TRUE)

或者像这样

x <- extract(r, p)
v <- sapply(x, mean)

每个多边形一个值

v
#[1] 14.00000 43.40000 49.00000 36.00000 29.50000 59.00000 91.00000 71.83333
#[9] 73.50000 87.16667 78.50000 59.57143

解决方案 2. 获得您要求的结构

x <- extract(r, p)
z <- do.call(rbind, lapply(1:length(x), function(i) cbind(i, x[[i]])))

或者像这样

i <- sapply(x, length)
j <- rep(1:length(i), i)
z <- cbind(j, unlist(x))

colnames(z) = c("ID", "value")
head(z)
#     ID value
#[1,]  1     4
#[2,]  1     5
#[3,]  1    13
#[4,]  1    14
#[5,]  1    15
#[6,]  1    22
tail(z)
#      ID value
#[51,] 12    55
#[52,] 12    56
#[53,] 12    57
#[54,] 12    64
#[55,] 12    65
#[56,] 12    66

稍后（见讨论中的问题）

要获取频率，您可以使用提取值列表和表函数。

使用新的示例值来获得频率的一些变化

values(r) <- rep(1:4, 25)

现在做

f <- extract(r, p, fun=function(i,...) table(i)) 

或分两步：

x <- extract(r, p) 
ff <- lapply(x, table)

每层都有一个表格

ff[1:2]
#[[1]]
#1 2 3 4 
#3 2 1 1 

#[[2]]
#1 2 3 4 
#1 1 2 1 

获得此结果的另一种方法是使用我上面所做的

x <- extract(r, p)
z <- do.call(rbind, lapply(1:length(x), function(i) cbind(i, x[[i]])))

紧随其后

fff <- tapply(z[,2], z[,1], table)

如果您想进一步操作这些，请为此编写一个函数并将其与lapply 或 for 循环一起使用。