如何将触摸多边形与 sf 包中的用户点数相等

How to unite touching polygons with equal number of users's points in sf package

我有一个带有大量非重叠多边形的 sf 对象,小样本是

polygons <- structure(list(id = c(138L, 150L, 138L, 138L, 153L), n_points = c(2L, 
3L, 3L, 2L, 1L), geometry = structure(list(structure(list(structure(c(2330, 
2380, 2380, 2371, 2330, 2330, 76, 76, 56, 56, 56, 76), .Dim = c(6L, 
2L))), class = c("XY", "POLYGON", "sfg")), structure(list(structure(c(2380, 
2380, 2371, 2371, 2380, 55, 46, 46, 55, 55), .Dim = c(5L, 2L))), class = c("XY", 
"POLYGON", "sfg")), structure(list(structure(c(2380, 2380, 2371, 
2371, 2380, 56, 55, 55, 56, 56), .Dim = c(5L, 2L))), class = c("XY", 
"POLYGON", "sfg")), structure(list(structure(c(2371, 2371, 2330, 
2330, 2371, 56, 55, 55, 56, 56), .Dim = c(5L, 2L))), class = c("XY", 
"POLYGON", "sfg")), structure(list(structure(c(2371, 2371, 2380, 
2380, 2330, 2330, 2371, 55, 46, 46, 40, 40, 55, 55), .Dim = c(7L, 
2L))), class = c("XY", "POLYGON", "sfg"))), n_empty = 0L, crs = structure(list(
    input = NA_character_, wkt = NA_character_), class = "crs"), class = c("sfc_POLYGON", 
"sfc"), precision = 0, bbox = structure(c(xmin = 2330, ymin = 40, 
xmax = 2380, ymax = 76), class = "bbox"))), row.names = c("138", 
"150", "138.1", "138.2", "153"), class = c("sf", "data.frame"
), sf_column = "geometry", agr = structure(c(id = NA_integer_, 
n_points = NA_integer_), .Label = c("constant", "aggregate", 
"identity"), class = "factor"))

> polygons
Simple feature collection with 5 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2330 ymin: 40 xmax: 2380 ymax: 76
CRS:           NA
       id n_points                       geometry
138   138        2 POLYGON ((2330 76, 2380 76,...
150   150        3 POLYGON ((2380 55, 2380 46,...
138.1 138        3 POLYGON ((2380 56, 2380 55,...
138.2 138        2 POLYGON ((2371 56, 2371 55,...
153   153        1 POLYGON ((2371 55, 2371 46,...

> plot(polygons["n_points"])

如何将点数相等的多边形连接起来?在我的示例中,它应该是 3 个多边形。 我可以使用 sf::st_touches(polygons) 找到接触的多边形,但我不知道下一步该做什么。

对于我的简单示例,期望的结果是

polygon_union <- split(polygons, polygons$n_points)
polygon_union <- purrr::map_dfr(polygon_union, dplyr::summarise)
polygon_union$n_points <- c(1, 2, 3)
plot(polygon_union["n_points"])

怎么样:

library(sf)
library(purrr)

num_points <- function (x) nrow(st_coordinates(x))

polygons$num_points <- map_dbl(seq_len(nrow(polygons)), 
                         ~num_points(polygons[.x,])
                       )

touches <- st_touches(polygons)

same_num_points <- map2(polygons$num_points, touches, 
                     function (x,y) x == polygons$num_points[y]
                   )

# filter touches
touches <- map2(touches, same_num_points, `[`)

# unite the touching polygons. 
# This will give warnings, as it throws away some feature information.

united <- map2(seq_len(nrow(polygons)), touches, 
      function(x, y) {
        x_polygon <- polygons[x,]
        y_polygons <- polygons[y,]
        st_union(x_polygon, y_polygons)
    })

# if you want these in a single data frame you can do:
united <- do.call(rbind, united)

请在下面找到一个仅使用 sf

的更简单的解决方案

Reprex

  • 代码
library(sf)

Results <- polygons %>% 
  aggregate(., 
            by = list(.$n_points), 
            function(x) x = floor(mean(x)), 
            join = st_touches)
  • 输出
Results
#> Simple feature collection with 3 features and 3 fields
#> Attribute-geometry relationship: 0 constant, 2 aggregate, 1 identity
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 2330 ymin: 40 xmax: 2380 ymax: 76
#> CRS:           NA
#>   Group.1  id n_points                       geometry
#> 1       1 153        1 POLYGON ((2371 55, 2371 46,...
#> 2       2 138        2 POLYGON ((2371 55, 2330 55,...
#> 3       3 144        3 POLYGON ((2380 46, 2371 46,...
  • 可视化
plot(Results["n_points"])

reprex package (v2.0.1)

于 2021-12-19 创建