R：制作加拿大邮政编码地图

我正在使用 R 编程语言。

这是我的问题的一些背景：

• 在加拿大，每个地址都有一个邮政编码（如美国邮政编码）。这些邮政编码的格式为：字母编号字母 - 数字字母编号（例如 A1B 2C3）
• 邮政编码的前三个字符称为 FSA - 所有前三个字母相同的邮政编码均被视为属于同一个 FSA
• 在线提供免费且公开的“地理空间形状文件”，其中包含加拿大所有 FSA 的地理边界（https://www150.statcan.gc.ca/n1/en/catalogue/92-179-X)
• 然而，加拿大邮政编码没有这样的形状文件 - 这就是我试图在 R 中构建的（即近似的）

为了开始解决这个问题，我找到了以下开放数据链接 - 其中包含数百万个经过地理编码的加拿大地址（每个地址的邮政编码和经度/纬度坐标）：https://www.statcan.gc.ca/en/lode/databases/oda

这是我的计划：

• 对于加拿大的某个省份 - 我可以创建一个“凸包”（https://en.wikipedia.org/wiki/Convex_hull) 具有相同邮政编码的所有地址。这可能可以作为每个邮政编码的“近似”边界。
• 如果我对所有邮政编码重复此操作，我将得到一个近似的形状文件

这是我到目前为止所尝试的：

Step 1:我下载了给定省份（例如安大略省）的文件（临时文件）：

url <- "https://www150.statcan.gc.ca/n1/en/pub/46-26-0001/2021001/ODA_ON_v1.zip"

temp <- tempfile()

download.file(url, temp)

unzipped_files <- unzip(temp, exdir = tempdir())

csv_file <- unzipped_files[grepl(".csv$", unzipped_files)]


df <- read.csv(unzipped_files[1], stringsAsFactors = FALSE)

Step 2:然后，使用lapply- 我使用内置函数计算了同一邮政编码中所有地址的凸包chull()功能。最后，我创建了一个形状文件，其最终结果是：

library(sf)

hulls <- lapply(unique(df$postal_code), function(postal_code) {
    chull(df[df$postal_code == postal_code, c("longitude", "latitude")])
})

hull_sfs <- lapply(seq_along(hulls), function(i) {
    st_as_sf(df[df$postal_code == unique(df$postal_code)[i], ][hulls[[i]], ],
             coords = c("longitude", "latitude"), crs = 4326)
})

hull_sf_combined <- do.call(rbind, hull_sfs)

st_write(hull_sf_combined, "hulls.shp")

Step 3:作为可选步骤，我们可以可视化最终结果：

library(ggplot2)

ggplot(df, aes(x = longitude, y = latitude)) +
    geom_point(aes(color = postal_code)) +
    lapply(seq_along(hulls), function(i) {
        geom_polygon(data = df[df$postal_code == unique(df$postal_code)[i], ][hulls[[i]], ],
                     aes(x = longitude, y = latitude, fill = postal_code), alpha = 0.5)
    })

我的问题：虽然代码似乎正在运行，但我不确定我所做的一切是否正确。有人可以帮我解决这个问题吗？

Thanks!

Note：根据@Ben Bolker提供的评论，是否最好采用此处提供的每个邮政编码的地理质心（https://www.serviceobjects.com/blog/free-zip-code-and-postal-code-database-with-geocooperatives/）并使用 Voronoi Cells 来确定近似边界？

这是有关清理邮政编码区域的部分答案。假设您使用凸包来获取一堆定义多边形的顶点。

library(polyclip)
library(ggplot2)

canada <- data.frame(x = c(0,1,1,0),
                     y = c(0,0,1,1))
postcode1 <- data.frame(x = c(0.5, 1.5, 1.5, 0.8),
                        y = c(0.5, 0.8, 1.5, 1.5))
postcode2 <- data.frame(x = c(0.25, 0.75, 0.75, 0.25),
                        y = c(0.25, 0.25, 0.75, 0.75))

ggplot() +
    geom_polygon(data = canada, aes(x,y), fill = "lightblue", color = "blue", alpha = 0.5) +
    geom_polygon(data = postcode1, aes(x,y), fill = "lightgreen", color = "green", alpha = 0.5) +
    geom_polygon(data = postcode2, aes(x,y), fill = "pink", color = "red", alpha = 0.5)

您可以看到邮政编码区域的凸包彼此重叠，其中之一延伸到国外。这polyclip包允许您对顶点所包围的区域进行交集、并集、异或和 setdiff。

clip1 <- polyclip(postcode1, canada, "intersect")[[1]] |> as.data.frame()
clip2 <- polyclip(postcode2, clip1, "minus")[[1]] |> as.data.frame()

ggplot() +
    geom_polygon(data = canada, aes(x,y), fill = "lightblue", color = "blue", alpha = 0.5) +
    geom_polygon(data = clip1, aes(x,y), fill = "lightgreen", color = "green", alpha = 0.5) +
    geom_polygon(data = clip2, aes(x,y), fill = "pink", color = "red", alpha = 0.5)

剪报快乐！