R GIS：基于邻居值添加多边形属性？

Question

我有一组分布随机的森林地块（SpatialPolygonDataFrame），即分散和聚集。对于每个多边形，我想判断它是否具有开放边缘。如果多边形具有开放边缘，则：

• 没有邻居
• 至少在一侧没有邻居;
• 有邻居，但是与邻居之间的树高度差大于5

我想知道如何将属性open_edge = TRUE / FALSE添加到各个多边形中？在raster包中，使用moving window可以有希望的方法。然而，我的原始数据是要素类，不像工作示例中那样类似于栅格。

我考虑过（伪代码）：

• 逐个子集化每个地块（在for循环中）
• 创建周围缓冲区
• 子集化周围地块与地块的buffer重叠部分
• 如果有任何邻居->比较高度。 如果差异> 5，则open_edge = TRUE

但是，这种方法没有考虑地块只有3个侧面邻居的情况，即车门邻域。 poly2nb工具似乎很有前途，但如何向各自地块添加属性呢？

这是我的虚拟方法，但我想知道您是否有更高效的解决方案？

创建虚拟数据：

library(ggplot2)  # for choropleth map plot
library(broom) # to convert spatial data to dataframe
library(mapproj)
library(spdep)    # neighbours
library(rgdal)
library(rgeos)
library(sf)
library(raster)
library(dplyr)
library(spData)
library(sf)

r <- raster(nrow=6, ncol=6, crs = "+init=epsg:2957")
values(r) <- matrix(data = c(9,  NA, NA, NA, NA, NA,
                             NA, NA, NA, NA, NA, NA, 
                             NA, NA, NA, NA, NA, NA, 
                             NA, NA, NA, 1, 1, 1,
                             NA, NA, NA, 1, 9, 1,
                             NA, NA, NA, 1, 1, 1),
                    nrow = 6,
                    ncol = 6, 
                    byrow = TRUE)

# Convert raster to polygon
polys <- rasterToPolygons(r)

通过一个展架示例来识别其是否具有开放的边缘：

# Subset first row in SpatialPolygonDataFrame
i = 10
one = polys[i, ]

# Keep the remaining polygons
left = polys[-i,]

# Create buffer within distance
buff = buffer(one, width = 100)

# subset set of neighbours by spatial overlap
nbrs <- left[which(gContains(sp::geometry(buff),
                                    sp::geometry(left), byid = T)),    
# Compare if the values are different
height.one  = rep(one$layer[1], nrow(nbrs))
height.nbrs = nbrs$layer

# Get the differences between the neighbouring stands
difference = height.one - height.nbrs

# If the difference in at least one stand is 
# in more than 5, set open_edge = TRUE 
# or if no neighbours find
one$open_edge <- any(difference > 5)

- maycca

Answer 1

让您开始使用spdep::poly2nb

library(raster)

r <- raster(nrow=6, ncol=6, crs = "+init=epsg:2957")
values(r) <- matrix(data = c(9,  NA, NA, NA, NA, NA,
                             NA, NA, NA, NA, NA, NA, 
                             NA, NA, NA, NA, NA, NA, 
                             NA, NA, NA, 1, 1, 1,
                             NA, NA, NA, 1, 9, 1,
                             NA, NA, NA, 1, 1, 1),
                    nrow = 6,
                    ncol = 6, 
                    byrow = TRUE)

# Convert raster to polygon
polys <- rasterToPolygons(r)

library(spdep)
nb <- poly2nb(polys)

plot(polys)
text(polys, 1:length(polys))

str(nb)
#List of 10
# $ : int 0
# $ : int [1:3] 3 5 6
# $ : int [1:5] 2 4 5 6 7
# $ : int [1:3] 3 6 7
# ...

因此，多边形1没有邻居，多边形2有邻居3、5、6等。

现在您可以使用sapply计算一些内容。例如，邻居的数量。

nbcnt <- sapply(nb, function(i) length(i[i>0]))
nbcnt 
#[1] 0 3 5 3 5 8 5 3 5 3

将此添加回多边形中。

polys$nbcnt <- nbcnt

- Robert Hijmans

---

感谢您的建议 @RobertHijmans。根据您的答案，我完成了循环遍历邻居以评估高度差异的函数。对于这种方法或如何简化它，您有什么想法？ - maycca

Answer 2

似乎有更简单的解决方案。您可以使用来自raster包的滑动窗口函数focal。

以下是一个示例：

library(raster)

r <- raster(nrow=6, ncol=6, crs = "+init=epsg:2957")
values(r) <- matrix(data = c(9,  NA, NA, NA, NA, NA,
                             NA, NA, NA, NA, NA, NA, 
                             NA, NA, 2, 1, 3, 1, 
                             NA, NA, 1, 1, 1, 1,
                             NA, NA, 1, 2, 2, 1,
                             NA, NA, 1, 1, 1, 1),
                    nrow = 6,
                    ncol = 6, 
                    byrow = TRUE)

# Prepare function for sliding window
focal_func <- function(x) {
    # Assuming 3x3 moving window
    # central cell has index 5
    # Check if the cell is not NA
    if (is.na(x[5])){
        return(NA)
    }

    # Check if any surrounding is NA
    if (any(is.na(x[-5]))) {
        return(TRUE)
    }

    # Check difference
    if (any((x[5] - x[-5]) > 5)) {
        return(TRUE)
    }

    # Else, return FALSE
    return(FALSE)
}

# Apply focal_function with sliding window
res <- focal(r, w = matrix(rep(1, 9), 3), fun = focal_func, pad = TRUE)

# Check if the same as desired output
res_mat <- as.matrix(res)
res[!is.na(res)] == 1

它给出：

[1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE  TRUE  TRUE FALSE FALSE
[13]  TRUE  TRUE  TRUE  TRUE  TRUE

即与所需输出相同。

- Istrel

---

谢谢您提供这种简洁明了的方法。不过，我将使用多边形而不是栅格，并且我不想进行转换。我的实际多边形可能具有奇怪的形状，类似于https://gis.stackexchange.com/questions/338679/spatial-subsetting-spatialpolygonsdataframes-completely-within-selecting-area-us - maycca

Answer 3

基于@RobertHijmans的答案和如何获取邻居列表，我创建了一个逐步检查我的邻居并评估高度差异的方法。

逐步操作：

1. 检查单元格有多少个邻居 如果邻居数量小于8，则将open_edge设置为TRUE，否则设置为FALSE
2. 循环遍历具有8个邻居（皇后情况）的单元格的索引： 比较中心和周围单元格之间的高度差异。 如果任何差异大于5，则将open_edge设置为TRUE。

这里是稍微复杂一些的情况，允许更多的支架有邻居：

创建虚拟数据：

r <- raster(nrow=6, ncol=6, crs = "+init=epsg:2957")
values(r) <- matrix(data = c(9,  NA, NA, NA, NA, NA,
                             NA, NA, NA, NA, NA, NA, 
                             NA, NA, 2, 1, 3, 1, 
                             NA, NA, 1, 1, 1, 1,
                             NA, NA, 1, 2, 2, 1,
                             NA, NA, 1, 1, 1, 1),
                    nrow = 6,
                    ncol = 6, 
                    byrow = TRUE)

# Convert raster to polygon
polys <- rasterToPolygons(r)

根据连续性计算邻居，考虑林分之间可能存在的间隙或细长部分：

# Calculate the count and position of neighbors; this allows to create one by one list
nb <- poly2nb(polys, 
              snap = 10) # snap corrects for the gaps/slivers

# store the number of neighbors for each central cell
polys$nb_count<- card(nb)

# Has the stand an open edge? Is surrounded by neighbors, pre-value is FALSE
polys$open_edge = ifelse(card(nb) <8, "TRUE", "FALSE")

如果一个单元格有完整的邻居，请比较它们之间高度的差异：

# Get the position of the cell surrounded by neighbors
center.index <- which(polys$nb_count == 8)

# Get the stand height of a stand
# as a vector to compare element wise
center.height = polys[center.index,]$layer

# Loop through the cells with neighbors:
# - keep height of the central stand
# - get height of neighbors
# compare the height between them
# if difference is more than 5: => open_edge = TRUE
for (i in seq_along(center.index)) {

  # Get central stand height 
  center.height = polys[center.index[i],]$layer

  # Identify neighbors of the stands
  # by the index value
  nb.index = unlist(nb[center.index[i]])

  # Get heights of the stands
  nb.height = polys[nb.index,]$layer

  # Adjust Center.height length as a vector to allow element wise comparison
  center.height.v = rep(center.height, length(nb.index))

  # Compare the heights 
  h.diff = center.height.v - nb.height

  # if any difference is more than 5 => change open_edge = TRUE
  if (any(h.diff > 5)) {
    polys@data[center.index[i],]$open_edge <- "TRUE"
  }
}

看最终的数据输出：

> polys@data
   layer nb_count open_edge
1      9        0      TRUE
2      2        3      TRUE
3      1        5      TRUE
4      3        5      TRUE
5      1        3      TRUE
6      1        5      TRUE
7      1        8     FALSE
8      1        8     FALSE
9      1        5      TRUE
10     1        5      TRUE
11     2        8     FALSE
12     2        8     FALSE
13     1        5      TRUE
14     1        3      TRUE
15     1        5      TRUE
16     1        5      TRUE
17     1        3      TRUE

- maycca