R光栅：：提取函数不按预期工作

Question

我试图在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)

但是，当多边形与栅格中的多个单元格相交时，这会返回多边形1的数字5。这是怎么回事？

然后我想得到加权平均值。如果有任何NA单元交叉，那么这应该‘失败’，并返回NA。我认为这是可行的，但由于上述错误，我现在对我的方法不太有信心。

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

，我也想知道相交的细胞数量少于一个，，我不知道该怎么做。有什么好主意吗？

polygons$almost_zero    <-

谢谢你的帮助。

Answer 1

来自extract帮助：

如果单元格的中心位于多边形内，则覆盖该单元格(但请参见考虑部分覆盖单元格的权重选项；对于获得小多边形值的参数小)。

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

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)

​

​

Answer 2

提取是按照它的意思，通过寻找细胞的质心在多边形内。处理部分覆盖的一种方法是使用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