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

Time:12-20

I have sf object with huge number of non-overlapping polygons, the small sample is

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"])

polygons

How to unite touching polygons with equal number of points? In my example its should be 3 polygons. I can find touching polygons using sf::st_touches(polygons) but I have no idea what to do next.

For my simple example the desired result is

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"])

result

CodePudding user response:

Please find below a more simple solution using only sf

Reprex

  • Code
library(sf)

Results <- polygons %>% 
  aggregate(., 
            by = list(.$n_points), 
            function(x) x = floor(mean(x)), 
            join = st_touches)
  • Output
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,...
  • Visualization
plot(Results["n_points"])

Created on 2021-12-19 by the reprex package (v2.0.1)

CodePudding user response:

How about:

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)
  •  Tags:  
  • r sf
  • Related