Leaflet和R中的投影坐标参考系统答案

【问题标题】：Projected coordinate reference systems in Leaflet and RLeaflet和R中的投影坐标参考系统
【发布时间】：2021-11-24 12:49:21
【问题描述】：

根据documentation，R 中的传单支持投影坐标参考系统。但是，我无法完成这项工作。

给定投影坐标系 EPSG:27700（称为 myprojectedsfobject）中的任何 sf 多边形对象，运行屏幕截图下方的代码块会给出代码中注释的错误。

crs <- leafletCRS(
                  code = "EPSG:27700",
                  proj4def =
                    "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs")

leaflet(options = leafletOptions(crs = crs)) %>%
  addPolygons(data = myprojectedsfobject, label = ~htmlEscape(NAME_3), fill = FALSE)
#Warning messages:
#1: sf layer is not long-lat data 
#2: sf layer has inconsistent datum (+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 #+y_0=-100000 +ellps=airy +units=m +no_defs).
#Need '+proj=longlat +datum=WGS84'

任何人都可以帮助如何在 R 的传单中正确实现投影坐标系吗？

【问题讨论】：

为什么不直接转换成长纬度数据呢？ myprojectedsfobject %>% st_transform(4326)
公平的挑战，但这是安全关键信息，其中由于转换导致的位置不准确可能是一个问题。另外，数据集很大，我不希望增加额外的处理开销。此外，我想与英国国家网格投影中的机载图像集成，这会导致另一个数据集的转换问题。简短的回答：设置投影坐标系会更容易。
我和 OP 有完全相同的问题，Jumble 的回答帮助解决了这个问题 - 已投赞成票。

标签： r leaflet gis

【解决方案1】：

从文档https://rstudio.github.io/leaflet/projections.html 看来，您仍然需要提供 WGS84 Lat Long 数据，然后传单将投影到您想要的 CRS。

来自上面的文档：

“虽然图块必须与 LeafletCRS 函数中使用的投影相同，但标记、圆、多边形和线必须始终使用 WGS 84 经度/纬度数据。Leaflet 将在显示时自动投影坐标。”

下面是 epsg:27700 瓦片和 epsg:4326 标记，它们会被传单自动转换为 27700。

epsg27700 <- leafletCRS(crsClass = "L.Proj.CRS", code = "EPSG:27700",
                       proj4def = "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs",
                       resolutions = c(896.0, 448.0, 224.0, 112.0, 56.0, 28.0, 14.0, 7.0, 3.5, 1.75),#2^(13:-1), # 8192 down to 0.5
                       origin = c(-238375.0, 1376256.0)
)


# EPSG:27700 TILE LAYER
tile_url <- "https://api.os.uk/maps/raster/v1/zxy/Leisure_27700/{z}/{x}/{y}.png?key=iJIgGUXKYFiT2s2ejG7cAB0sGuvtOCyp"

leaflet(options = leafletOptions( crs = epsg27700)) %>%
  addTiles(urlTemplate = tile_url)%>%
  
  # WGS84 Markers are transformed to EPSG:27700 behind the scenes
  addMarkers(-0.16359, 51.5083)

【讨论】：