在R中基于sf包为不同区域的坐标点集绘制外围多边形

需求背景 #

有多个来自不同地理区域的点簇数据集，需要为每个区域生成单个外围多边形，用于地图绘制空间范围，对几何计算精度要求低。

示例测试数据 #

以下代码可生成结构一致的测试数据，包含挪威、美国本土、阿拉斯加三个区域的点数据：

# 测试数据生成代码
nor_lat <- c(72.117, 71.05, 71.717, 71.167, 71.817, 72.333, 71.25, 70.917, 
  70.85, 70.933)
nor_lon <- c(33.217, 30.333, 26.3, 22.333, 22.333, 22.333, 18.517, 20.133, 
             21.333, 19.917)
us_lat <- c(41.283333, 41, 40.95, 40.95, 38.683333, 40.783333, 40.733333, 
            40.516667, 38.566667, 41.266667)
us_lon <- c(-71.116667, -70.733333, -70.583333, -70.483333, -74.816667, 
            -70.5, -70.566667, -70.266667, -74.85, -71.366667)
ak_lat <- c(58.015, 57.67167, 57.33833, 56.685, 56.99333, 57.31667, 57.65, 
            57.99667, 58.32167, 58.33333)
ak_lon <- c(-158.31667, -158.36, -158.41, -159.76, -159.72333, -159.66333, 
            -159.64167, -159.605, -159.54, -160.72167)
nor_dat <- data.frame(lat=nor_lat, lon=nor_lon)
us_dat <- data.frame(lat=us_lat, lon=us_lon)
ak_dat <- data.frame(lat=ak_lat, lon=ak_lon)
nor_dat$region <- "nor"
us_dat$region <- "us"
ak_dat$region <- "ak"
dat <- rbind(nor_dat, us_dat, ak_dat)

实现步骤 #

依赖包加载 #

需要用到sf做空间处理，dplyr做分组批量操作：

library(sf)
library(dplyr)

转换为空间点对象 #

将普通数据框转换为sf格式的空间点对象，指定WGS84坐标参考系（CRS:4326）：

dat_sf <- st_as_sf(dat, coords = c("lon", "lat"), crs = 4326)

单组点生成外围多边形 #

以挪威区域为例，单组点生成外围轮廓的代码如下：

# 筛选挪威区域的点集
nor_points <- filter(dat_sf, region == "nor")
# 合并点集后生成凸包（即外围轮廓多边形）
nor_poly <- st_convex_hull(st_union(nor_points))

执行完成后直接运行plot(nor_poly)即可查看生成的多边形效果。

批量生成所有区域的外围多边形 #

无需循环，直接通过分组聚合即可一次性生成所有区域的外围多边形：

region_polys <- dat_sf %>%
  group_by(region) %>%
  summarise(geometry = st_union(geometry)) %>%
  st_convex_hull()

生成的region_polys是包含三个区域多边形的sf对象，可直接用于绘图，还可以叠加原始点验证轮廓正确性：

# 绘制分区域的外围多边形
plot(region_polys["region"], key.pos = 1)
# 叠加原始点
plot(st_geometry(dat_sf), pch = 16, col = "red", add = TRUE)

可选优化 #

如果需要更贴合点分布的平滑凹轮廓，可以替换st_convex_hull()为concaveman包的concaveman()函数，仅用于绘图的场景下默认的凸包方案已足够使用。

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