R raster::extract 函数未按预期工作答案

【问题标题】：R raster::extract function not working as expectedR raster::extract 函数未按预期工作
【发布时间】：2018-05-21 10:54:44
【问题描述】：

我正在尝试在 R 中做一些多边形/栅格工作。提取功能没有像我预期的那样工作。这是一个可重现的示例。

library(sf)
library(raster)
library(maptools)

raster      <-     download.file('https://raw.githubusercontent.com/JimShady/london_osm_canyons/master/test.tif', 'downloaded_raster')
raster      <- raster('test.tif')

polygons    <- read_sf('https://raw.githubusercontent.com/JimShady/london_osm_canyons/master/test2.geojson')
polygons    <- st_transform(polygons, crs(raster)@projargs)
polygons$id <- 1
polygons    <- as(polygons, 'Spatial')

new_polygon <- elide(polygons, rotate=-50)
proj4string(new_polygon) <- crs(raster)@projargs

polygons    <- rbind(polygons, new_polygon)

plot(raster)
plot(polygons, add=T)

我想获取每个多边形相交的单元格的数量，即使它们是 NA，并将它们作为列存储在多边形中

polygons$cell_count     <- extract(raster, polygons, fun = function(x,...)length(x), na.rm=F)

但是，当多边形明显与栅格的 5 个以上像元相交时，这将返回第一个多边形的数字 5。 这是怎么回事？

然后我想得到加权平均值。如果有任何 NA 单元格相交，那么这应该“失败”并返回 NA。 我认为这行得通，但是由于上面的错误，我现在对我的方法不是很有信心。

polygons$mean           <- round(extract(raster, polygons, fun = mean, weights=T, na.rm=F),2)

我还想知道相交小于 1 的单元格的数量。我不知道该怎么做。有人有什么好主意吗？

polygons$almost_zero    <-

感谢您的帮助。

【问题讨论】：

运行您的上述代码时出现许多错误。确保代码易于重现。
对不起，我认为它有效。我现在就修。
现在应该可以工作了@SamAct

标签： r gis raster r-raster sf

【解决方案1】：

来自extract帮助：

如果一个单元格的中心在多边形内，则该单元格被覆盖（但请参见考虑部分覆盖单元的权重选项；和论据 small 无论如何获取小多边形的值）。

因此，您可以使用它来获得单元格的数量：

data <- extract(raster, polygons,  weights = T,  na.rm = F)[[1]]
nrow(data)

#[1] 18

data 在哪里：

      value      weight
 [1,]     8 0.003610108
 [2,]     6 0.025270758
 [3,]     0 0.010830325
 [4,]     7 0.104693141
 [5,]     0 0.066787004
 [6,]     0 0.088447653
 [7,]     0 0.001805054
 [8,]     0 0.104693141
 [9,]     0 0.010830325
[10,]     0 0.093862816
[11,]     0 0.030685921
[12,]     0 0.063176895
[13,]     0 0.057761733
[14,]     0 0.063176895
[15,]     7 0.003610108
[16,]     7 0.102888087
[17,]     0 0.093862816
[18,]     0 0.074007220

拥有相交且小于 1 的单元格数：

sum(data[, 1] < 1)

#[1] 13

编辑

如果有多个多边形，您应该首先使用rgeos::gUnaryUnion 分解您的多边形：

library(rgeos)

polygons$id <- 1
polygons <- gUnaryUnion(polygons, id = polygons$id)
# As suggested by @RobertHijmans, you can also use `raster::aggregate`
# polygons <- raster::aggregate(polygons, by = "id")

plot(raster)
plot(polygons, add = T)

【讨论】：

谢谢 DJack。但是，如果有多个多边形，我将如何处理？我已经编辑了我的问题以反映这种情况。我想也许我需要使用某种 lapply 但我正在努力解决它。谢谢。
@TheRealJimShady 您不应该以使现有答案不再正确的方式真正编辑问题。如果您有进一步的后续问题，您应该将它们作为新问题发布（或者如果它们是微不足道的，或者只是对当前答案的澄清）
@TheRealJimShady 我进行了编辑。不过dww是对的，下次如果有后续问题，应该发个新问题。
欢呼并注意到
要“溶解”，使用aggregate 更方便，因为它保留了属性

【解决方案2】：

通过查找质心位于多边形内的单元格，Extract 的行为符合预期。处理部分覆盖的一种方法是使用rasterize 将多边形转换为光栅蒙版。请注意，在下面我分别重命名了您的栅格和多边形 r1 和 p1。因为对对象使用函数名并不是一个好习惯：

r_mask = rasterize(p1, r1, getCover=TRUE)
r_mask[r_mask==0] = NA

现在我们可以使用这个掩码来获取你想要的值：

r2 = mask(r1, r_mask)
cellStats(!is.na(r2), sum)
# [1] 18
cellStats(r2, mean)
# [1] 1.944444

要找出有多少个单元格小于一，我们可以这样做

cellStats(r2<1, sum)
# [1] 13

【讨论】：