如何在R中将等高线合并为多边形时保留孔洞？

你遇到的问题是st_polygonize将等高线转换为多边形时，无法自动识别内环（孔洞），导致中心区域被填充。以下是两种在R中解决该问题的方法：

方法1：手动调整等高线环的方向构造带孔洞多边形 #

rasterToContour生成的等高线包含两条闭合线（外环和内环），但st_polygonize会将它们视为独立多边形。我们可以通过调整环的方向（外环与内环方向相反），手动构造带孔洞的多边形：

library(raster)
library(sf)

# 生成示例栅格
r <- raster(matrix(c(rep(0,5),
                     0,1,1,1,0,
                     0,1,0,1,0,
                     0,1,1,1,0,
                     rep(0,5)), nrow = 5, ncol = 5))

# 生成等高线并转为sf对象
contour_lines <- rasterToContour(r, level = 0.5)
sf_lines <- st_as_sf(contour_lines)

# 提取每条线的坐标
line_coords <- lapply(sf_lines$geometry, function(x) st_coordinates(x)[,1:2])

# 判断环的方向（面积正负区分顺/逆时针）
get_orientation <- function(coords) {
  area <- sum((coords[-1,1] - coords[-nrow(coords),1]) * (coords[-1,2] + coords[-nrow(coords),2])) / 2
  ifelse(area > 0, "counter-clockwise", "clockwise")
}
orientations <- sapply(line_coords, get_orientation)

# 确定外环（面积更大的环）
outer_idx <- which.max(sapply(line_coords, function(x) abs(sum((x[-1,1]-x[-nrow(x),1])*(x[-1,2]+x[-nrow(x),2]))/2)))
outer_coords <- line_coords[[outer_idx]]
inner_coords <- line_coords[[which(orientations != orientations[outer_idx])]]

# 调整外环为逆时针，内环为顺时针
if (orientations[outer_idx] == "clockwise") outer_coords <- outer_coords[nrow(outer_coords):1,]
if (get_orientation(inner_coords) == "counter-clockwise") inner_coords <- inner_coords[nrow(inner_coords):1,]

# 构造带孔洞的多边形
poly_with_hole <- st_polygon(list(outer_coords, inner_coords))
sf_poly <- st_sfc(poly_with_hole, crs = st_crs(sf_lines)) %>% st_sf()

# 验证绘图
plot(r)
plot(sf_poly, col = "gray50", border = "black", add = TRUE)

方法2：直接从栅格生成带孔洞多边形（更简便） #

如果你的需求是基于栅格值生成带孔洞的区域多边形，直接使用raster::rasterToPolygons更高效，它会自动识别内部的孔洞区域：

library(raster)
library(sf)

r <- raster(matrix(c(rep(0,5),
                     0,1,1,1,0,
                     0,1,0,1,0,
                     0,1,1,1,0,
                     rep(0,5)), nrow = 5, ncol = 5))

# 从栅格生成多边形，dissolve=TRUE合并相同值的区域
poly_raster <- rasterToPolygons(r, dissolve = TRUE)
sf_poly <- st_as_sf(poly_raster)

# 筛选出值为1的带孔洞多边形
sf_poly_1 <- sf_poly[sf_poly$layer == 1,]

# 验证绘图
plot(r)
plot(sf_poly_1, col = "gray50", border = "black", add = TRUE)

第二种方法无需手动处理等高线，适合直接从栅格值提取区域的场景。

内容的提问来源于stack exchange，提问作者DDaye