如何用R语言sf包将坐标生成简单轮廓多边形而非乱线？

问题描述 #

现有一组无序坐标，尝试用sf包直接转换为多边形时，生成的图形是点间杂乱连线而非包含所有点的轮廓。原始代码如下：

library(sf)
library(tidyverse)
lon <- c(-71.616667, -71.5, -71.716667, -71.9, -72.4, -72.65, -72.383333, 
-72.4, -72.383333, -72.65, -72.6, -72.7, -73.033333, -73.2, -73.933333, 
-73.9, -74, -74.1, -74.183333, -74.183333, -74.45, -74.566667, 
-74.533333, -73.55, -73.733333, -73.833333, -73.683333, -73.283333, 
-72.883333, -73, -69.966667, -69.8, -73.183333, -73.85, -73.533333, 
-73.983333, -74.1, -74.3, -73.916667, -73.983333, -73.933333, 
-74.183333, -74.35, -74.366667, -74.516667, -74.583333, -73.283333, 
-73.15, -72.8, -72.583333)
lat <- c(39.916667, 39.916667, 39.85, 39.95, 39.816667, 39.766667, 39.633333, 
39.6, 39.316667, 39.266667, 39.15, 39.116667, 39.05, 39.033333, 
39.983333, 39.816667, 39.7, 39.516667, 39.35, 39.283333, 39.116667, 
39.066667, 39.033333, 39.2, 39.333333, 39.333333, 39.616667, 
39.616667, 39.566667, 39.85, 39.95, 39.933333, 39.8, 39.683333, 
39.333333, 39.15, 39.15, 39.066667, 39.85, 39.733333, 39.65, 
39.266667, 39.216667, 39.183333, 39.1, 39.033333, 39.1, 39.15, 
39.45, 39.383333)

tmp <- data.frame(lon, lat)

# 原始转换代码（生成杂乱图形）
polygon <- tmp %>%
    st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
    summarise(geometry = st_combine(geometry)) %>%
    st_cast("POLYGON")

问题根源 #

原始坐标是无序的，st_cast("POLYGON")会直接按点的原始顺序连线，导致图形混乱。需要先计算包含所有点的轮廓边界，最简便的方式是生成凸包（Convex Hull）——即包含所有点的最小凸多边形。

解决方案 #

1. 生成凸包多边形 #

修改转换代码，用st_convex_hull()替代直接st_cast("POLYGON")：

# 生成包含所有点的凸包多边形
polygon <- tmp %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_convex_hull()  # 关键：自动计算凸包边界

# 查看结果
plot(polygon)

2. 精确计算面积（可选） #

因为原始坐标是WGS84（EPSG:4326，地理坐标系），直接计算的面积是球面近似值，建议转换为投影坐标系（如UTM）后再计算：

# 转换为北美东部UTM投影（EPSG:32618，根据坐标位置选择合适分区）
polygon_utm <- polygon %>% st_transform(crs = 32618)

# 计算面积（单位：平方米）
area_m2 <- st_area(polygon_utm)
print(area_m2)

3. 非凸边界需求（可选） #

如果需要更贴合的非凸边界，可先对点进行顺时针/逆时针排序，再生成多边形：

• 计算点集的中心坐标
• 按点与中心的角度排序
• 再转换为多边形

示例代码：

# 计算点集中心
center <- tmp %>% summarise(lon = mean(lon), lat = mean(lat))

# 计算每个点与中心的角度并排序
tmp_sorted <- tmp %>%
  mutate(angle = atan2(lat - center$lat, lon - center$lon)) %>%
  arrange(angle)

# 转换为多边形
polygon_sorted <- tmp_sorted %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON")

plot(polygon_sorted)

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