区域统计 R

Posted 2023-02-23

【中文标题】区域统计 R

【英文标题】：Zonal Statistics R

【发布时间】：2016-06-03 13:40:16

【问题描述】：

我在R 中使用raster 包并在shapefile 包含多个polygons (200) 的光栅上执行extract。每个多边形都有一个列，其中指定了多边形的class。 当我执行extract 时，我得到一个dataframe，其中栅格的每个值都“分配”给一个多边形：哪个多边形位于哪个像素中。但是，我的多边形包含的信息（即它们是什么类）已经消失了。 extract 之后的我的数据框看起来像：

test <- extract(myRaster, myTrainingPolygon, df = TRUE)

ID     band1     band2    band3
1      0.101     0.827    ...
...    ...       ...      ...
200    0.876     0.821    ...

我需要的是类似的东西

ID     band1     band2    band3    class
1      0.101     0.827    ...      class1
...    ...       ...      ...      ...
200    0.876     ...      ...      class3

我如何在其中获取信息 - 或者首先 - 在执行提取时不要丢失它们？！

【问题讨论】：

你能提供一个你的 shapefile 的例子吗？

【参考方案1】：

始终包含并使用示例数据

library(raster)
r <- raster(ncol=36, nrow=18)
r[] <- 1:ncell(r)
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))
p <- spPolygons(cds1, cds2)
p <- SpatialPolygonsDataFrame(p, data.frame(class=c('A', 'B'), stringsAsFactors=FALSE))

现在我们可以这样做了：

v <- extract(r, p)
str(v)

#List of 2
# $ : num [1:38] 326 327 328 329 330 331 332 333 334 335 ...
# $ : num [1:25] 172 173 208 209 244 245 279 280 281 282 ...

创建一个具有多边形顺序 ID 和提取值的 data.frame

d <- data.frame(id=rep(1:length(v), sapply(v, length)), value=unlist(v))

（这种使用 unlist 的特殊方法仅适用于单层）。对于多层做

#d <- data.frame(id=rep(1:length(v), sapply(v, length)), do.call(rbind, v))

使用顺序 ID 和所需的其他多边形属性创建一个 data.frame

pd <- cbind(id=1:length(p), data.frame(p))

合并两者

m <- merge(pd, d)

【参考方案2】：

?extract 表示sp 参数控制是否应将提取的值添加到空间对象的数据框中。指定 sp=TRUE 应该可以解决问题。例如：

library(maptools)
library(raster)
data(wrld_simpl)
r <- raster(extent(-180, 180, -90, 90), res=10)
r[] <- runif(ncell(r))
wrld_simpl_new <- extract(r, wrld_simpl, fun=mean, sp=TRUE)

head(wrld_simpl_new)

##     FIPS ISO2 ISO3 UN                NAME   AREA  POP2005 REGION SUBREGION     LON     LAT     layer
## ATG   AC   AG  ATG 28 Antigua and Barbuda     44    83039     19        29 -61.783  17.078 0.9142067
## DZA   AG   DZ  DZA 12             Algeria 238174 32854159      2        15   2.632  28.163 0.1774097
## AZE   AJ   AZ  AZE 31          Azerbaijan   8260  8352021    142       145  47.395  40.430 0.3098710
## ALB   AL   AL  ALB  8             Albania   2740  3153731    150        39  20.068  41.143 0.3746480
## ARM   AM   AM  ARM 51             Armenia   2820  3017661    142       145  44.563  40.534 0.3494729
## AGO   AO   AO  AGO 24              Angola 124670 16095214      2        17  17.544 -12.296 0.2873931

添加的列名为“layer”，因为这是RasterLayer 的（默认）名称。

【讨论】：

这需要将fun 参数设置为某些东西，例如mean。但是，我只想知道哪些像素位于哪个栅格中。假设您有 100x100 像素的栅格，并且有两个多边形。在一个多边形中有 50 个像素，在其他 20 个中。我想写出哪些像素在哪个多边形中！。执行mean 会将我的所有像素减少到一个值。