如何基于点数据近似加拿大邮政编码区多边形并生成Shapefile

数据加载（已完成） #

你已通过以下代码加载了包含地址、经纬度和邮政编码的数据集：

library(httr)
library(data.table)

url <- "https://www.dropbox.com/scl/fi/9kjoqsppb85ip0tdc5wmr/stackoverflow_example.csv?rlkey=cwjk222jnoz8c9cgbt01p6bep&dl=1"
response <- GET(url)
df <- fread(content(response, "text"), sep = ",", quote = "", fill = TRUE)

需求目标 #

利用现有地址点数据，近似每个邮政编码的地理边界，最终导出包含边界经纬度的Shapefile，同时接受其他可行方案。

方案一：多边形裁剪（Polygon Clipping）实现边界生成 #

该方法通过泰森多边形（Voronoi）+ 行政边界裁剪生成贴合实际的邮政编码边界，是你提到的Polygon Clipping技术的典型应用：

步骤与代码 #

1. 加载空间处理依赖包

library(sf)
library(dplyr)
library(tmap)

2. 将数据转换为空间点对象（采用WGS84坐标系，EPSG:4326）

sf_points <- st_as_sf(df, coords = c("longitude", "latitude"), crs = 4326)

3. 生成泰森多边形（初步划分点集边界）

# 生成全局Voronoi多边形并转换为可处理格式
voronoi_polys <- st_voronoi(st_union(sf_points)) %>% 
  st_cast("POLYGON") %>% 
  st_sf()

# 将Voronoi多边形与原始地址点关联，匹配对应邮政编码
voronoi_postal <- voronoi_polys %>% 
  st_join(sf_points, join = st_contains) %>% 
  group_by(postal_code) %>% 
  summarise(geometry = st_union(geometry))  # 合并同一邮编的多边形

4. 用行政边界裁剪Voronoi多边形（避免边界超出实际城市范围）

# 从现有数据中提取Kingston的近似行政边界（实际可替换为官方边界数据）
city_boundary <- sf_points %>% 
  filter(csdname == "Kingston") %>% 
  st_union() %>% 
  st_convex_hull()

# 执行多边形裁剪
clipped_postal <- voronoi_postal %>% 
  st_intersection(city_boundary)

5. 可视化验证与导出Shapefile

# 可视化边界与原始点
tm_shape(clipped_postal) +
  tm_polygons(col = "postal_code", alpha = 0.5) +
  tm_shape(sf_points) +
  tm_dots(size = 0.1, col = "black")

# 导出为Shapefile
st_write(clipped_postal, "postal_code_clipped_boundaries.shp", delete_layer = TRUE)

方案二：凸包（Convex Hull）快速近似边界 #

如果不需要高精度边界，凸包是更轻量化的方案，直接生成每个邮编点集的最小外包多边形：

# 按邮政编码分组生成凸包
postal_hulls <- sf_points %>% 
  group_by(postal_code) %>% 
  summarise(geometry = st_convex_hull(st_union(geometry)))

# 可视化
tm_shape(postal_hulls) +
  tm_polygons(col = "postal_code", alpha = 0.5) +
  tm_shape(sf_points) +
  tm_dots(size = 0.1)

# 导出Shapefile
st_write(postal_hulls, "postal_code_hull_boundaries.shp", delete_layer = TRUE)

方案三：官方边界数据（最优精度） #

若需要精确的邮政编码边界，建议直接使用加拿大统计局的官方数据，可通过cancensus包获取（需申请API密钥）：

install.packages("cancensus")
library(cancensus)

# 示例：获取指定区域的邮政编码边界（需替换为你的目标区域）
postal_official <- get_census(
  dataset = "CA21",
  regions = list(CSD = "3510010"),  # Kingston的CSD代码
  level = "postal",
  vectors = "C1_001",
  api_key = "你的API密钥"  # 从cancensus官网申请
)

# 导出官方边界Shapefile
st_write(postal_official, "postal_code_official_boundaries.shp", delete_layer = TRUE)

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