在 R 中使用矢量多边形提取栅格像素值
Extract Raster Pixels Values Using Vector Polygons in R
我已经为此苦苦挣扎了几个小时。
我有一个包含 177 个多边形(即 177 个县)的 shapefile(称为 "shp")。此 shapefile 覆盖在栅格上。我的栅格(称为 "ras")由具有不同污染值的像素组成。
现在我想提取每个多边形的所有像素值及其出现次数。
这正是 QGIS 功能 "zonal histogram" 所做的。但我想在 R 中做完全相同的事情。
我尝试了 extract() 函数,我设法得到了每个县的平均值,这已经是第一步,但我想制作像素分布(直方图)。
有人可以帮帮我吗?
非常感谢,
玛丽·劳尔
这是一个最小的、独立的、可重现的例子(几乎字面上来自 ?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]))
非常感谢您的帮助。下次我保证我会小心并更详细地解释我的问题。
在您的帮助下,我设法找到了解决方案。
我也用过这个网站: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)
我已经为此苦苦挣扎了几个小时。 我有一个包含 177 个多边形(即 177 个县)的 shapefile(称为 "shp")。此 shapefile 覆盖在栅格上。我的栅格(称为 "ras")由具有不同污染值的像素组成。
现在我想提取每个多边形的所有像素值及其出现次数。
这正是 QGIS 功能 "zonal histogram" 所做的。但我想在 R 中做完全相同的事情。
我尝试了 extract() 函数,我设法得到了每个县的平均值,这已经是第一步,但我想制作像素分布(直方图)。
有人可以帮帮我吗?
非常感谢,
玛丽·劳尔
这是一个最小的、独立的、可重现的例子(几乎字面上来自 ?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]))
非常感谢您的帮助。下次我保证我会小心并更详细地解释我的问题。
在您的帮助下,我设法找到了解决方案。 我也用过这个网站: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)