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