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)

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
# ...

所以 poly 1 没有邻居，poly 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

Answer 2

似乎有一个更简单的解决方案。您可以使用 focal 包中的 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

即与所需的输出相同。

Answer 3

根据@RobertHijmans 的回答和how to get list of neighbors，我创建了一个逐步检查邻居和评估身高差异的方法。

一步一步：

1. 检查单元有多少邻居 如果邻居数 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