在R中为分组经纬度坐标生成外边界多边形的技术问询

你的原代码直接将所有坐标点按输入顺序连接成多边形，这不是真正的"最外侧点外边界"，会出现交叉、凹陷的不规则形状，完全不符合空间分析需求。要生成精准的外边界，推荐用凸包（Convex Hull）（适合规则形状设施）或阿尔法形状（Alpha Shape）（适合非凸、有凹陷的复杂设施），以下是适配你需求的完整实现：

1. 加载依赖包 #

library(sf)
library(dplyr)
library(alphashape3d) # 用于生成阿尔法形状，需先安装：install.packages("alphashape3d")

2. 读取并预处理数据 #

假设你的数据已加载为df，先转换为sf点对象：

# 转换为sf点图层（替换为你数据中的经纬度列名）
sf_points <- st_as_sf(df, coords = c("longitude", "latitude"), crs = 4326) # WGS84坐标系

3. 自定义分组处理函数 #

该函数会根据每组建筑数量，返回单个点或精准外边界多边形：

generate_boundary <- function(geom) {
  # 过滤空几何
  geom_filtered <- geom[!st_is_empty(geom)]
  n_points <- length(geom_filtered)
  
  if (n_points == 0) {
    return(st_sfc(st_point(), crs = 4326)) # 返回空点
  } else if (n_points == 1) {
    return(geom_filtered) # 保留单个建筑的点几何
  } else {
    # 优先生成阿尔法形状（贴合点分布的非凸边界）
    coords <- st_coordinates(geom_filtered)
    # alpha参数需按需调整：值越小边界越贴合点，越大越接近凸包
    alpha_opt <- ashape(coords, alpha = 0.01)$alpha[which.max(ashape(coords, alpha = 0.01)$area)]
    shape <- ashape(coords, alpha = alpha_opt)
    
    # 若阿尔法形状生成失败（如点过于密集）， fallback到凸包
    if (length(shape$edges) == 0) {
      hull <- st_convex_hull(st_union(geom_filtered))
      return(hull)
    } else {
      # 将阿尔法形状转换为sf多边形
      edges <- shape$edges
      edge_coords <- coords[c(edges$ind1, edges$ind2), ]
      polygon <- st_polygon(list(rbind(edge_coords, edge_coords[1, ]))) %>% st_sfc(crs = 4326)
      return(polygon)
    }
  }
}

4. 分组执行并生成结果 #

result <- sf_points %>%
  group_by(installation_name) %>%
  summarise(geometry = generate_boundary(geometry)) %>%
  ungroup()

5. 导出为GIS兼容格式 #

直接导出为可导入ArcGIS/QGIS的文件：

# 导出GeoJSON（推荐，无文件大小限制）
st_write(result, "facility_boundaries.geojson", driver = "GeoJSON")

# 或导出Shapefile
st_write(result, "facility_boundaries.shp", driver = "ESRI Shapefile")

关键说明 #

• 凸包vs阿尔法形状：凸包是所有点的最小外接凸多边形，计算速度快但会忽略设施内部凹陷；阿尔法形状能生成贴合建筑分布的非凸边界，适合有内部空地的大型设施，alpha参数可根据数据测试调整（比如0.005、0.01、0.05）。
• 效率优化：针对数万条数据，sf的st_union和凸包计算效率极高；若遇到50+建筑的大型设施，可适当调大alpha值提升阿尔法形状的生成速度。
• 坐标系统：示例使用WGS84（EPSG:4326），若需更高精度的空间分析，可通过st_transform(result, crs = 326XX)转换为UTM平面坐标系（XX替换为对应带号）。

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