无法在tmap与leaflet中绘制OSM多面几何体问题排查

问题描述 #

我用osmdata从OpenStreetMap提取数据，查询得到包含点、线、多边形及多面几何体的数据集。使用plot(st_geometry(...))可正常绘制该多面几何体，但tmap和leaflet却无法识别并绘制它们，请问这是什么原因？

复现代码 #

library(osmdata)
library(sf)
library(tmap)

nkpg_schools <- 
  opq(bbox = "Norrköping") %>% 
  add_osm_feature(key = "amenity", value = "school") %>%
  osmdata_sf()

plot(st_geometry(nkpg_schools$osm_multipolygons)) # 正常显示

tmap_mode("view")
tm_basemap("Esri.WorldGrayCanvas") +
  tm_shape(nkpg_schools$osm_multipolygons) +
  tm_polygons()
# 无法显示

原因分析与解决方法 #

核心原因 #

1. 坐标系不兼容：tmap的交互式视图默认采用Web墨卡托投影（EPSG:3857），而osmdata返回的sf对象默认是WGS84地理坐标系（EPSG:4326）。虽然tmap理论上支持自动转换，但复杂MultiPolygon可能因转换失败导致渲染异常。
2. 几何拓扑错误：OSM中的部分MultiPolygon可能存在自相交、无效环等拓扑问题，plot()函数对这类错误的容忍度较高，但tmap/leaflet依赖的底层渲染库对几何有效性要求更严格。
3. 空几何对象干扰：提取的结果中可能包含空的MultiPolygon，会导致渲染流程中断。

分步解决 #

1. 手动转换坐标系

在传入tmap前，显式将sf对象转换为Web墨卡托投影：

# 转换坐标系
nkpg_schools_proj <- st_transform(nkpg_schools$osm_multipolygons, crs = 3857)

# 重新绘图
tmap_mode("view")
tm_basemap("Esri.WorldGrayCanvas") +
  tm_shape(nkpg_schools_proj) +
  tm_polygons()

2. 修复几何拓扑错误

使用sf的st_make_valid()函数修复无效的几何对象：

# 修复几何有效性
nkpg_schools_valid <- st_make_valid(nkpg_schools$osm_multipolygons)
# 转换坐标系
nkpg_schools_valid_proj <- st_transform(nkpg_schools_valid, crs = 3857)

# 绘图
tmap_mode("view")
tm_basemap("Esri.WorldGrayCanvas") +
  tm_shape(nkpg_schools_valid_proj) +
  tm_polygons()

3. 过滤空几何对象

先检查并移除结果中的空几何：

# 过滤空几何
nkpg_schools_filtered <- nkpg_schools$osm_multipolygons[!st_is_empty(nkpg_schools$osm_multipolygons), ]
# 转换坐标系后绘图
nkpg_schools_filtered_proj <- st_transform(nkpg_schools_filtered, crs = 3857)

tmap_mode("view")
tm_basemap("Esri.WorldGrayCanvas") +
  tm_shape(nkpg_schools_filtered_proj) +
  tm_polygons()

内容的提问来源于stack exchange，提问作者Maël