我有一个地球上的矢量多边形列表,有些是独立且不同的,有些是接触的。通常,在现实中,感人的部分都是同一形状的一部分,但以这种方式分割,因为不同的年份添加了不同的部分,等等。
我需要以编程方式组合那些正在接触的,但并非全部都是,有时许多连续接触,端到端。
我认为应该有一种方法可以使用
relate()
和 combineGeoms()
或类似的方法来做到这一点,但除了手动编辑每个形状组之外,我不知道如何做到这一点。对于大型且复杂的形状文件来说,这似乎是一个可行的任务。
也欢迎使用 sf 的解决方案,我不介意哪种方式。
我想要的是一个函数,我可以使用关系矩阵来合并组内接触的所有形状,从而产生行数较少的矢量数据帧,因为单独的接触形状将分解为一个。在这种情况下,数据集将从 7 个多边形变为 4 个多边形。
下面的最小代表:
library(sf)
library(dplyr)
library(ggplot2)
# Define a function to create a simple square polygon
create_square <- function(x_center, y_center, size) {
half_size <- size / 2
st_polygon(list(matrix(c(
x_center - half_size, y_center - half_size,
x_center + half_size, y_center - half_size,
x_center + half_size, y_center + half_size,
x_center - half_size, y_center + half_size,
x_center - half_size, y_center - half_size
), ncol = 2, byrow = TRUE)))
}
# Create individual shapes
shape1 <- create_square(1, 1, 1)
shape2 <- create_square(4, 1, 1)
shape3 <- create_square(6, 1, 1)
shape4 <- create_square(7, 1, 1)
shape5 <- create_square(7.5, 2, 1)
shape6 <- create_square(2.5, 2.5,1)
shape7 <- create_square(2, 3.5, 1)
# Combine shapes into an sf object
shapes <- st_sf(
geometry = st_sfc(
shape1,
shape2,
shape3,
shape4,
shape5,
shape6,
shape7
)
) %>%
mutate(name = paste0("shape_",c(1:nrow(.)))) %>%
relocate(name)
# Define a coordinate reference system (CRS), for example, EPSG:4326 (WGS84)
crs <- st_crs(4326)
st_crs(shapes) <- crs
# plot
shapes %>%
ggplot() +
geom_sf(aes(fill = name))
# convert to spatvector
library(terra)
library(tidyterra)
shapes_vect <- vect(shapes)
ggplot() +
geom_spatvector(data = shapes_vect,
aes(fill = name))
# make relation matrix of touching shapes
relate(shapes_vect,
shapes_vect,
relation = "touches")
# manually combine shapes 3,4,5
test <- combineGeoms(shapes_vect[3,],
shapes_vect[5,]) %>%
combineGeoms(shapes_vect[4,])
# create new dataframe with the combined shape
# instead of 3,4,5 individually
shapes_new <- shapes_vect[c(1:2, 6:7),] %>%
rbind(test)
# plot again
ggplot() +
geom_spatvector(data = shapes_new,
aes(fill = name))
dplyr::summarise()
与 sf::st_union()
结合将达到您想要的结果。请注意,将一些 sf
函数应用于非投影数据,例如WGS84/EPSG:4326,可能会导致不需要的伪影。
为了解决这个问题,您可以首先投影数据,然后在多边形合并后将其更改回原始 CRS。在此示例中,我使用 WGS84 Pseudo Mercator/EPSG:3857 来实现此目的。
library(sf)
library(dplyr)
library(ggplot2)
# Project shapes, merge contiguous polygons, st_cast() to individual polygons,
# change back to original CRS
sf_poly <- shapes |>
st_transform(3857) |>
summarise(geometry = st_union(geometry)) |>
st_cast("POLYGON") |>
st_transform(4326)
sf_poly
# Simple feature collection with 4 features and 0 fields
# Geometry type: POLYGON
# Dimension: XY
# Bounding box: xmin: 0.5 ymin: 0.5 xmax: 8 ymax: 4
# Geodetic CRS: WGS 84
# geometry
# 1 POLYGON ((1.5 1.5, 1.5 0.5,...
# 2 POLYGON ((4.5 0.5, 3.5 0.5,...
# 3 POLYGON ((5.5 0.5, 5.5 1.5,...
# 4 POLYGON ((2 2, 2 3, 1.5 3, ...