绘制 ggmap 和 sf 点

Posted 2023-03-23

【中文标题】绘制 ggmap 和 sf 点

【英文标题】：plot ggmap and sf points

【发布时间】：2022-01-17 02:46:46

【问题描述】：

我正在努力从 ggmap 绘制带有点（sf 对象）的底图。我遇到的最有希望的解决方案是来自@andyteucher 的SO answer。我试图在下面重现它，但我运气不佳。我得到了要打印的底图，但没有得到点。

library(tidyverse)
library(sf)
library(ggmap)

# Define a function to fix the bbox to be in EPSG:3857
# https://***.com/a/50844502/841405
  ggmap_bbox <- function(map) 
    if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
    # Extract the bounding box (in lat/lon) from the ggmap to a numeric vector, 
    # and set the names to what sf::st_bbox expects:
    map_bbox <- setNames(unlist(attr(map, "bb")), 
                         c("ymin", "xmin", "ymax", "xmax"))
    
    # Coonvert the bbox to an sf polygon, transform it to 3857, 
    # and convert back to a bbox (convoluted, but it works)
    bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), 3857))
    
    # Overwrite the bbox of the ggmap object with the transformed coordinates 
    attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
    attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
    attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
    attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
    map
  
# requires API key
  basemap <- get_map(location=c(lon = 75.85199398072335, 
                                lat = 22.7176905515565), 
                     zoom=9, maptype = 'toner-hybrid', 
                     source = 'google')


  basemap_3857 <- ggmap_bbox(basemap)

  points <- tribble(
    ~name, ~lat, ~lon,
    "test1", 22.7176905515565, 75.85199398072335,
    "test2", 22.71802612842761, 75.84848927237663,
  ) %>%
    st_as_sf(coords = c("lat", "lon"),
             crs = 3857)

  ggmap(basemap_3857) + 
    coord_sf(crs = st_crs(3857)) + 
    geom_sf(data = points, 
            inherit.aes = FALSE) 

【问题讨论】：

你好@Eric Green。我不使用ggmap，所以很难帮助你。也就是说，在我看来，您的代码中有一个小错误。在我看来coords = c("lat", "lon") 应该是coords = c("lon", "lat")。希望这可以帮助。干杯。

谢谢@lovalery。不幸的是，它并不能解决我最初的问题。

这只是一个没有太多信念的尝试 ;-) 不幸的是我不能更好地帮助你，因为我不是 ggmap 的专家。祝你好运！

【参考方案1】：

我相信您的坐标参考系统有问题 - 您似乎在 CRS 3857 的上下文中使用了度数，它以米为单位定义（因此相差几个度数......）

如果我的预感是正确的，您需要首先将您的 sf 对象声明为 4326（WGS84 = GPS 坐标的 CRS），然后 - 并且只有在那时 - 将转换应用到 3857（从已知开始）。

此代码和地图符合您的期望吗？ （我还稍微清理了 get_map 调用，因为它混合了 Google 和 Stamen 术语；没什么大不了的）

# requires API key
basemap <- get_map(location=c(lon = 75.85199398072335, 
                              lat = 22.7176905515565), 
                   zoom=9, 
                   source = 'google')


basemap_3857 <- ggmap_bbox(basemap)

points <- tribble(
  ~name, ~lat, ~lon,
  "test1", 22.7176905515565, 75.85199398072335,
  "test2", 22.71802612842761, 75.84848927237663,
) %>%
  st_as_sf(coords = c("lon", "lat"),
           crs = 4326) %>% # this is important! first declare wgs84
  st_transform(3857) # and then transform to web mercator 

ggmap(basemap_3857) + 
  geom_sf(data = points,
          color = "red",
          inherit.aes = F) 

【讨论】：

直到时间的尽头，我对 crs 和转换的看法都是错误的

@EricGreen :D 这是一个容易出错的话题 - 我与亚丁湾建立了深厚而亲密的关系，这是我通常感兴趣的领域（捷克***）出现的地方当旋转 90° 时（= 经度交换）。大约一半的时间发生在我身上，即完全随机