R中Shapefile导入后Leaflet可视化报错及解决原理咨询

问题描述 #

我尝试用R的sf库导入Shapefile：

library(sf)
library(leaflet)

sf <- sf::st_read("C:/Users/me/OneDrive/Documents/shape5/shapefile.shp", options = "ENCODING=WINDOWS-1252") 

导入输出信息：

options:        ENCODING=WINDOWS-1252 
Reading layer `lda_000a16a_e' from data source `C:\Users\me\OneDrive\Documents\shape5\shapefile.shp' using driver `ESRI Shapefile'
Simple feature collection with 56590 features and 22 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 3658201 ymin: 658873 xmax: 9019157 ymax: 6083005
Projected CRS: PCS_Lambert_Conformal_Conic

用Leaflet绘制时出现警告且无可视化结果：

> leaflet()  %>% addTiles() %>%  addPolygons(data=sf,weight=5,col = 'red')

警告信息：

Warning messages:
1: sf layer is not long-lat data 
2: sf layer has inconsistent datum (+proj=lcc +lat_0=63.390675 +lon_0=-91.8666666666667 +lat_1=49 +lat_2=77 +x_0=6200000 +y_0=3000000 +datum=NAD83 +units=m +no_defs).
Need '+proj=longlat +datum=WGS84' 

后来用以下代码处理后可视化正常：

sf <- sf::st_read("C:/Users/me/OneDrive/Documents/shape5/shapefile.shp", options = "ENCODING=WINDOWS-1252") %>%    
  sf::st_transform('+proj=longlat +datum=WGS84')

原因解析 #

• 坐标系类型不匹配：Leaflet地图默认采用WGS84经纬度坐标系（即+proj=longlat +datum=WGS84），用经度（-180_{180）和纬度（-90}90）表示位置。而你的Shapefile用的是Lambert Conformal Conic投影坐标系，坐标单位是米，属于平面投影坐标，Leaflet无法直接解析这种非经纬度格式的坐标，因此报“not long-lat data”警告，无法渲染图形。
• 地理基准面不一致：你的数据采用NAD83基准面，和Leaflet要求的WGS84基准面是两种不同的地理参考框架，二者对地球形状的建模存在细微差异，坐标数值对应的实际地理位置会有偏差，所以提示“inconsistent datum”。
• st_transform的核心作用：sf::st_transform()函数的功能就是空间坐标系转换，它会将Shapefile数据从原来的Lambert投影+NAD83基准面，精准转换为Leaflet要求的WGS84经纬度坐标系。转换后的数据坐标格式、基准面完全匹配Leaflet的要求，自然就能正常显示可视化效果。

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