我有一堆代表土地覆盖类别的多边形,但是我不需要这么多细节。我想合并小多边形(即< 150m2) with it's largest neighbour. This would be similar to "Eliminate" in ArcGIS Pro or "eliminate selected polygons" in QGIS. This process will be repeated so I'd like to script it in R and save some time.
我能想到的唯一方法是添加一个列来指示哪些多边形小于 150m2,然后以某种方式将它们合并到最近的邻居。
#Packages
library(terra)
library(sf)
library(dplyr)
#Adding polygons from Terra and converting to sf
v <- vect(system.file("ex/lux.shp", package="terra"))
v <- st_as_sf(v[c(1:12)])
#Adding a column indicating which polygons are under 150m2
mutate(v, area_thresh = case_when(AREA < 150 ~ 1,
AREA > 150 ~ 0))
这是一个有趣的问题,因为您需要遍历不断变化的 shapefile 行(随着合并较小的对象,您的行数会减少)。
对于一种可能的方法,请考虑这段代码,该代码建立在众所周知且备受喜爱的 NC shapefile 上,随
{sf}
.
请注意,我使用 county_id 作为唯一标识符;合并的多边形将保留较大县的 ID。
library(sf)
library(dplyr)
# the one & only NC shapefile... as geometry only
shape <- st_read(system.file("shape/nc.shp", package="sf")) %>%
select(CNTY_ID)
# limit for merging - in this case median area
limit <- median(st_area(shape))
# initialize iterator
i <- 1
# iterate over rows of shapefile
while (i < nrow(shape)) {
# check if area over or under limit
if(st_area(shape[i, ]) < limit) {
# find id of largest neighbor
largest_neighbor <- shape[unlist(st_touches(shape[i, ], shape)), ] %>%
mutate(area = st_area(.)) %>%
top_n(1, area) %>%
pull(CNTY_ID)
# merge offending shape with largest neighbor
wrk_shape <- st_union(st_geometry(shape[i, ]),
st_geometry(shape[shape$CNTY_ID == largest_neighbor, ])) %>%
st_make_valid() %>%
st_as_sf() %>%
mutate(CNTY_ID = largest_neighbor)
# rename geometry, column to enable binding
st_geometry(wrk_shape) = attr(shape, "sf_column")
# remove offending shape and unmerged neighbor
shape <- shape[-i, ]
shape <- filter(shape, !CNTY_ID == largest_neighbor)
# add merged shape back in
shape <- bind_rows(shape, wrk_shape)
}
# don't forget to update iterator :)
i <- i + 1
}
plot(st_geometry(shape))