在R/Leaflet重现LTE覆盖地图时遇孤立孔洞错误的解决咨询

目标 #

使用R语言在Leaflet中重现2019年加拿大LTE覆盖地图。

已执行步骤 #

步骤1：下载JSON数据 #

library(jsonlite)
library(httr)

response <- GET("https://crtc.gc.ca/cartovista/LTEOverTheYearsYE2019_EN/map/LTE_YE2019.json")
content <- content(response, as = "text")
data <- fromJSON(content)

步骤2：提取多边形并设置CRS #

library(sp)
library(rgdal)

polygons <- lapply(data$f$g$c[[1]], function(x) {
  Polygon(matrix(x, ncol = 2, byrow = TRUE))
})

sp_polygons <- SpatialPolygons(list(Polygons(polygons, ID = 1)))
spdf <- SpatialPolygonsDataFrame(sp_polygons, data.frame(id = 1, row.names = 1))
proj4string(spdf) <- CRS(data$proj)

数据结构（截取部分）：

> str(spdf)
Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 1 obs. of  1 variable:
  .. ..$ id: num 1
  ..@ polygons   :List of 1
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 606
  ....
  .. .. .. ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. .. ..@ projargs: chr "+proj=lcc +lat_0=49 +lon_0=-95 +lat_1=49 +lat_2=77 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"

遇到的错误 #

使用Leaflet绘图时出现拓扑错误：

library(leaflet)

m <- leaflet() %>%
  addTiles() %>%
  addPolygons(data = spdf)

# Error in rgeos::createPolygonsComment(pgons) : 
#   rgeos_PolyCreateComment: orphaned hole, cannot find containing polygon for hole at index 2

尝试过的修复方法 #

方法1：使用cleangeo清理空间对象 #

require(devtools)
install_github("eblondel/cleangeo")
require(cleangeo)

report <- clgeo_CollectionReport(spdf)
summary <- clgeo_SummaryReport(report)
issues <- report[report$valid == FALSE,]

spdf.clean <- clgeo_Clean(spdf)

方法2：使用rmapshaper简化多边形 #

library(rmapshaper)

spdf_clean <- ms_simplify(spdf)

m <- leaflet() %>%
  addTiles() %>%
  addPolygons(data = spdf_clean)

方法3：使用sf包修复拓扑有效性 #

library(sf)
library(leaflet)

sf_polygons <- st_as_sf(spdf)
valid <- st_is_valid(sf_polygons)
sf_polygons_valid <- st_make_valid(sf_polygons)

m <- leaflet() %>%
  addTiles() %>%
  addPolygons(data = sf_polygons_valid)

解决方案 #

1. 处理路径正确性评估 #

你的基础处理路径是正确的：获取原始JSON数据→解析为空间对象→设置正确CRS。错误核心是多边形存在孤立孔洞，这是空间数据常见拓扑问题，你的三种修复方向均合理，但可优化执行效率与准确性。

2. 高效解决孤立孔洞的优化方案 #

优先使用sf包工作流，它的拓扑修复工具更现代高效，且Leaflet对sf对象原生支持。优化步骤如下：

步骤1：直接用sf解析并处理数据

跳过sp包手动构建多边形的步骤，减少中间错误：

library(sf)
library(httr)
library(jsonlite)
library(leaflet)

# 获取并解析数据
response <- GET("https://crtc.gc.ca/cartovista/LTEOverTheYearsYE2019_EN/map/LTE_YE2019.json")
content <- content(response, as = "text")
data <- fromJSON(content)

# 转换为sf多边形对象
polygon_list <- lapply(data$f$g$c[[1]], function(coords) {
  st_polygon(list(matrix(coords, ncol = 2, byrow = TRUE)))
})

sf_obj <- st_sfc(polygon_list, crs = data$proj)
sf_df <- st_sf(id = 1:length(polygon_list), geometry = sf_obj)

# 修复拓扑问题（自动处理孤立孔洞）
sf_valid <- st_make_valid(sf_df)

# 转换为Leaflet支持的WGS84坐标系（EPSG:4326）
sf_valid_4326 <- st_transform(sf_valid, 4326)

步骤2：绘制Leaflet地图

m <- leaflet() %>%
  addTiles() %>%
  addPolygons(data = sf_valid_4326, 
              fillColor = "#1E90FF", 
              fillOpacity = 0.3, 
              color = "#1E90FF", 
              weight = 1)
m

3. 方案优势 #

• 减少手动错误：直接用sf的st_polygon构建对象，避免sp包中孔洞与父多边形的匹配错误。
• 针对性拓扑修复：st_make_valid自动将孤立孔洞转换为独立多边形，比其他工具的修复更精准。
• Leaflet原生兼容：sf对象无需额外转换即可被识别，且补充了坐标转换步骤（原始代码缺少此关键环节）。

4. 额外优化建议 #

• 简化数据提升加载速度：修复后可使用rmapshaper简化多边形顶点数量：

library(rmapshaper)
sf_simplified <- ms_simplify(sf_valid_4326, keep = 0.1, keep_shapes = TRUE)

• 匹配原地图样式：调整addPolygons的配色参数，还原原地图的蓝色覆盖层样式。

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