在 ggplot2 中绘制底图（basemaps_gglayer() 不起作用）

Question

我正在尝试为我的一些图形添加底图。从网上阅读来看，似乎

basemaps

包是最好的方法，但我一直无法让它工作，并且遇到了几个问题。

这是该图的原始代码，运行得很好：

ggplot() +
  geom_raster(subset(results, !is.na(var)), mapping = aes(long, lat, fill = var), interpolate = TRUE) +
  coord_sf(datum = "ESRI:102001") +
  scale_fill_gradientn('Variance', colours = colors) + #limits = c(0, 1)) +
  geom_sf(data = lg.parks, fill = NA, color = "black") + 
  geom_sf(data = canada, fill = NA, color = "black") + 
  annotation_scale(location = "bl", width_hint = 0.5) +
  annotation_north_arrow(location = "bl", which_north = "true", 
                         pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
                         style = north_arrow_fancy_orienteering) +
  theme_void() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank())

这是我现在的身材：

我想为此图添加底图。

当我尝试添加行

basemap_gglayer(ext = ext(canada)) + 

时，它返回错误

Error in st_transform.sfc(st_as_sfc(st_bbox(x)), crs = crs_webmerc) : 
  cannot transform sfc object with missing crs

我似乎可以通过指定

ext = st_bbox(canada)

来解决这个问题。但是，如果我尝试用

basemap_gglayer(ext = st_bbox(canada))

线绘制我的图形，我会收到另一个错误：

Error in `scale_fill_gradientn()`:
! Discrete values supplied to continuous scale.
ℹ Example values: "#91DED8", "#91DED8", "#91DED8", "#91DED8", and "#91DED8"

所以这是行不通的。我尝试将底图保存为单独的栅格，然后使用

geom_spatraster

包中的

tidyterra

函数来绘制它。

bm <- basemaps::basemap_raster(ext = st_bbox(canada))
bm <- as(bm, "SpatRaster")

ggplot() +
  geom_spatraster(data = bm) +
  geom_raster(subset(results, !is.na(var)), mapping = aes(long, lat, fill = var), interpolate = TRUE) +
  coord_sf(datum = "ESRI:102001") +
  scale_fill_gradientn('Variance', colours = colors) + #limits = c(0, 1)) +
  geom_sf(data = lg.parks, fill = NA, color = "black") + 
  geom_sf(data = canada, fill = NA, color = "black") + 
  annotation_scale(location = "bl", width_hint = 0.5) +
  annotation_north_arrow(location = "bl", which_north = "true", 
                         pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
                         style = north_arrow_fancy_orienteering) +
  theme_void() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank())

这是输出：

这真的不是我想要的，哈哈。如果我只绘制基础层，它看起来像这样：

理想情况下，我的最终图形看起来像^这样，我的原始图形只是分层在上面。

我的默认值设置为

basemaps::set_defaults(map_service = "esri", map_type = "world_terrain_base")

，因为我想避免使用需要 API 密钥的 google/stadia 地图。我正在 linux ubuntu v 22.04.4 上运行 R 4.3.3。

我在下面生成了一些可重现的代码：

library(ggplot2)
library(rgeoboundaries) #for canada shapefile
library(basemaps)
library(tidyterra)
library(terra) #for crop + mask functions

#canada shapefile
canada <- geoboundaries("Canada")

basemaps::set_defaults(map_service = "esri", map_type = "world_terrain_base")

#generate a raster
x <- raster(ncol=50, nrow=50, xmn=-141, xmx=-52, ymn=41, ymx=83)
y <- crop(x, canada)
values(y) <- 1:ncell(y)
y <- mask(y, canada)
z <- as.data.frame(y, xy = TRUE) #convert to dataframe for ggplot

ggplot() +
  #basemap_gglayer(ext = st_bbox(canada)) + 
geom_raster(subset(z, !is.na(layer)), mapping = aes(x, y, fill = layer), interpolate = TRUE) +
  coord_sf(datum = "ESRI:102001") +
  geom_sf(data = canada, fill = NA, color = "black") + 
  theme_void()

我不确定这是包裹的问题还是我的问题。或者也许有更好的方法来获取底图！任何建议表示赞赏

Answer 1

您的方法存在一些问题：

• ESRI

  basemap_raster()

  数据的默认 CRS 为 EPSG:3857，但

  geoboundaries()

  数据的 CRS 为 EPSG:4326；
• 对于具有 RGB 值的 SpatRaster，您正在使用

  geom_spatraster()

  而不是

  geom_spatraster_rgb()

  ；
• 您的 x 对象没有 CRS，这对于空间对象来说通常不是一个好的做法，并且 x 没有扩展到您的 canada 对象的完整范围（小问题，但会导致地图难看;)）。

以下代表通过为所有数据设置一致的 CRS、扩展 x 的范围以及使用正确的绘图函数来解决所有这些问题。请注意，您没有提供 colors 和 lg.parks 对象的副本，因此这些已被替换或省略。

加载包并创建数据：

library(rgeoboundaries) # For Canada shapefile
library(basemaps)
library(tidyterra)
library(sf) # For projecting sf objects
library(terra)
library(ggplot2)
library(ggspatial)

# Canada shapefile, project to match default basemap_raster crs
canada <- geoboundaries("Canada") %>%
  st_transform(3857)

# Basemap
set_defaults(map_service = "esri", map_type = "world_terrain_base")
bm <- basemap_raster(ext = st_bbox(canada)) %>%
  as("SpatRaster")

# Get extent of Canada sf for min/max values of example x raster
st_bbox(canada)
#      xmin      ymin      xmax      ymax 
# -15696344   5113338  -5857561  17955343 

# Generate an example raster using canada min/max values
x <- rast(ncols = 50, nrows = 50, nlyrs = 1, 
          xmin = -15696344, xmax = -5857561, 
          ymin = 5113338, ymax = 17955343, 
          names = c("rast_val"), 
          crs = "EPSG:3857") %>%
  init("cell")

y <- crop(x, canada)
values(y) <- 1:ncell(y)
y <- mask(y, canada)
z <- as.data.frame(y, xy = TRUE) # Convert to dataframe for ggplot

绘图数据：

ggplot() +
  geom_spatraster_rgb(data = bm) +
  geom_raster(subset(z, !is.na(rast_val)),
              mapping = aes(x, y, fill = rast_val, alpha = 0.5),
              interpolate = TRUE) +
  geom_sf(data = canada, fill = NA, color = "black") +
  scale_fill_gradientn("Variance", colours = terrain.colors(5)) +
  # geom_sf(data = lg.parks, fill = NA, color = "black") + 
  annotation_scale(location = "bl", 
                   width_hint = 0.5,
                   pad_x = unit(0.18, "in")) +
  annotation_north_arrow(location = "bl",
                         which_north = "true", 
                         pad_x = unit(0.2, "in"), 
                         pad_y = unit(4, "in"),
                         style = north_arrow_fancy_orienteering) +
  guides(alpha = "none") +
  coord_sf(datum = "ESRI:102001") +
  theme_void()

将会出现以下警告：

地图上的比例尺误差超过10%，比例尺可能不准确

结果：