设置 CRS 以使用雄蕊图进行绘图

我试图在地图上绘制一个简单的多边形来表示我感兴趣的区域。迄今为止，我已将多边形定义为并且能够单独绘制它。

poly <- st_polygon(list(as.matrix(data.frame(lat = c(40, 40, 60, 60, 40),
                                             lon = c(-60, -40, -40, -60, -60))))) #EDIT: these are lat/lon coordinates
ggplot() + 
    geom_sf(data = poly,
    aes(geometry = geometry),
    col = "red)

但是，我需要添加底图，以便可以显示该多边形相对于其他空间元素的位置。

# Attempt #1: stamenmaps basemap

    grid_box <- c(left   = -70, 
                  right  = -30,
                  bottom = 35,
                  top    = 70)
    base <- get_stamenmap(grid_box, zoom = 6, maptype = "terrain-background", color = "color")
    ggmap(base) +
      geom_sf(data = poly, 
              aes(geometry = geometry),
              col = "red")

    Coordinate system already present. Adding new coordinate system, which will replace the existing one.
    Error in `geom_sf()`:
    ! Problem while computing aesthetics.
    ℹ Error occurred in the 4th layer.
    Caused by error in `FUN()`:
    ! object 'lon' not found
    Run `rlang::last_error()` to see where the error occurred.

我发现向多边形添加 crs 的唯一方法如下（st_transform 和 st_as_sf 不起作用），但这极大地改变了坐标的比例，并且仍然无法绘制。

 # Attempt #2: new CRS

    poly <- st_polygon(list(as.matrix(data.frame(lon = c(-62, -43, -43, -62, -62),
                                                     lat = c(43, 43, 70, 70, 43))))) %>%
            st_sfc(crs = 3857)
    ggmap(base) +
      geom_sf(data = poly, 
              aes(geometry = geometry),
              col = "red")

Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Error in `geom_sf()`:
! Problem while computing aesthetics.
ℹ Error occurred in the 4th layer.
Caused by error in `FUN()`:
! object 'lon' not found
Run `rlang::last_error()` to see where the error occurred.

How can I get this polygon to plot over my basemaps?

我相信你走在正确的道路上sf::st_sfc()请再次尝试。为了让它发挥作用，你必须想一下你的坐标意味着什么？

In 每股收益：3857 https://epsg.io/3857坐标以米为单位定义。这似乎不是您名为 lat & long 的原始列的情况，它们通常代表纬度和经度（以度为单位）。

如果您的坐标以度为单位（这似乎很可能，尽管我在您的问题中找不到明确说明），那么使用 EPSG:4326 会更好。

您必须解决的一个技巧是 ggmap 被标记为 EPSG:4326（未投影度），但实际上是在 EPSG:3857（墨卡托投影米）中绘制的，这里描述了一种可能的方法如何将 geom_sf 生成的地图放在 ggmap 生成的栅格之上 https://stackoverflow.com/questions/47749078/how-to-put-a-geom-sf-produced-map-on-top-of-a-ggmap-produced-raster

基于所链接的问题，我提出以下建议：

library(sf)
library(ggmap)

poly <- st_polygon(list(as.matrix(data.frame(lon = c(-62, -43, -43, -62, -62),
                                             lat = c(43, 43, 70, 70, 43))))) %>%
  st_sfc(crs = 4326) 

grid_box <- c(left   = -70, 
              right  = -30,
              bottom = 35,
              top    = 70)

base <- get_stamenmap(grid_box, zoom = 5, 
                      maptype = "terrain-background", 
                      color = "color")

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
}

# Use the function:
base <- ggmap_bbox(base)

ggmap(base) +
  coord_sf(crs = st_crs(3857)) +
  geom_sf(data = st_transform(poly, 3857), 
          col = "red", fill = NA,
          inherit.aes = F) +
  theme(axis.title = element_blank())