在 R 中使用矢量多边形提取光栅像素值答案

【问题标题】：Extract Raster Pixels Values Using Vector Polygons in R在 R 中使用矢量多边形提取光栅像素值
【发布时间】：2020-04-02 13:30:18
【问题描述】：

我已经为此苦苦挣扎了好几个小时。我有一个包含 177 个多边形（即 177 个县）的 shapefile（称为“shp”）。此 shapefile 覆盖在栅格上。我的栅格（称为“ras”）由具有不同污染值的像素组成。

现在我想提取每个多边形的所有像素值及其出现次数。

这正是 QGIS 功能“区域直方图”正在做的事情。但我想在 R 中做同样的事情。

我尝试了 extract() 函数，并设法获得了每个县的平均值，这已经是第一步，但我想制作一个像素分布（直方图）。

有人可以帮帮我吗？

非常感谢，

玛丽-劳尔

【问题讨论】：

标签： r gis raster

【解决方案1】：

非常感谢您的帮助。下次我保证我会小心并更详细地解释我的问题。

在您的帮助下，我设法找到了解决方案。我也用过这个网站：http://zevross.com/blog/2015/03/30/map-and-analyze-raster-data-in-r/

关于信息，首先我必须卸载“tidyr”包，因为与提取功能有冲突。

如果它可以帮助某人，这是最终代码：

# Libraries loading
library(raster) 
library(rgdal)
library(sp)

# raster layer import
ras=raster("C:/*.tif")

# shapefile layer import
shp<-shapefile("C:/*.shp")

# Extract the values of the pixels raster per county
ext <- extract(ras, shp, method='simple')

# Function to tabulate pixel values by region & return a data frame
tabFunc                            <- function(indx, extracted, region, regname) {
  dat                              <- as.data.frame(table(extracted[[indx]]))
  dat$name                         <- region[[regname]][[indx]]
  return(dat)
}

# run through each county & compute a table of the number
# of raster cells by pixel value. ("CODE" is the county code) 
tabs <- lapply(seq(ext), tabFunc, ext, shp, "CODE")

# assemble into one data frame
df <- do.call(rbind, tabs)  

# to see the data frame in R
print(df)

# table export 
write.csv(df,"C:/*.csv", row.names = FALSE)

【讨论】：

【解决方案2】：

这是一个最小的、独立的、可重现的示例（几乎完全来自?raster::extract，所以不难制作）

library(raster)
r <- raster(ncol=36, nrow=18, vals=rep(1:9, 72))
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- spPolygons(cds1, cds2)

现在你可以做

v <- extract(r, polys)
par(mfrow=c(1,2))
z <- lapply(v, hist)

或者更花哨

mains <- c("first", "second")
par(mfrow=c(1,2))
z <- lapply(1:length(v), function(i) hist(v[[i]], main=mains[i]))

或者你想要一个条形图

z <- lapply(1:length(v), function(i) barplot(table(v[[i]]), main=mains[i]))

【讨论】：