如何获取sf对象的有效数据外包络轮廓？

我正在寻找一种从sf对象中提取/推断出多边形的方法，该多边形需描述包围有效数据的轮廓。我使用sf的经验有限，因此除了解答外，任何相关建议也都很受用。

示例数据（特殊投影） #

假设我有一些采用特殊投影的数据：

library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE

world <- st_as_sf(rworldmap::getMap(resolution = "low"))

tworld <- st_transform(world, "+proj=igh")
plot(tworld, max.plot = 1)

我可以顺利获取其格网（不过对于IGH投影而言，生成的格网可能不够美观）：

grat <- st_graticule(tworld)
plot(grat, max.plot = 1)

但如何可靠地获取包围数据的轮廓呢？我曾考虑直接重新投影一个经纬度矩形，但当初始几何图形并非完整世界地图时，这种方法是否依然有效？如果初始数据只是笛卡尔XY坐标，或许我应该直接获取X和Y的范围来构建矩形：

p <- st_polygon(
  list(matrix(c(-180, -180, 180, 180, -180, -90, 90, 90, -90, -90), 5))
)
p <- st_segmentize(p, 1)
p <- st_sfc(p)
st_crs(p) <- 4326 # 经纬度坐标
proj_p <- st_transform(p, "+proj=igh")
plot(proj_p)

解决方案与相关建议 #

方法1：轴对齐最小包围矩形（适合笛卡尔投影） #

如果数据是平面笛卡尔投影，直接用边界框生成多边形是最简单高效的方式：

# 获取sf对象的边界框
bbox <- st_bbox(tworld)
# 将边界框转换为sf多边形
bbox_poly <- st_as_sfc(bbox)
# 可视化验证
plot(tworld, max.plot = 1)
plot(bbox_poly, border = "red", lwd = 2, add = TRUE)

注意：对于IGH这类分瓣投影，轴对齐矩形会覆盖大量空白区域，仅适合数据集中在单一区域的场景。

方法2：凸包（Convex Hull） #

凸包能生成包含所有要素的最小凸多边形，比轴对齐矩形更贴合数据整体轮廓：

# 合并所有要素为单个几何对象
combined <- st_union(tworld)
# 计算凸包
convex_hull <- st_convex_hull(combined)
# 可视化
plot(tworld, max.plot = 1)
plot(convex_hull, border = "blue", lwd = 2, add = TRUE)

局限：凸包会忽略数据中的凹进部分，如果你的数据有明显的凹陷区域，结果会不够精准。

方法3：阿尔法形状（Alpha Shape） #

阿尔法形状可以生成更贴合数据细节的多边形，甚至保留凹进结构，适合复杂轮廓的场景：
首先需要安装alphashape3d包：

install.packages("alphashape3d")
library(alphashape3d)

# 提取所有要素的顶点坐标
points <- st_coordinates(st_cast(tworld, "POINT"))
# 计算阿尔法形状，alpha值需根据数据调整（越小越贴合细节）
alpha_shape <- ashape(points, alpha = 0.1)
# 转换为sf格式的多边形
alpha_sf <- st_as_sfc(alpha_shape, crs = st_crs(tworld))
# 可视化
plot(tworld, max.plot = 1)
plot(alpha_sf, border = "green", lwd = 2, add = TRUE)

关键提示：alpha值需要反复调试，不同数据的最优值差异较大。

特殊投影（如IGH）的处理注意事项 #

• 全局经纬度矩形投影到分瓣投影（如IGH）会得到不连续的多边形，仅适用于完整世界地图，局部数据绝对不能用这种方法。
• 处理局部数据时，优先使用边界框、凸包或阿尔法形状，避免跨投影的全局范围转换。
• 所有操作必须保证对象处于同一投影坐标系下，否则会出现几何计算错误。

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