使用terra绘图与maptiles地图时add_grid函数异常问题

我通过maptiles库配合terra的SpatRaster创建带背景地图的可视化，初始代码可以生成正常的网格：

library(terra)
library(maptiles)

# 示例栅格：局部高程数据
r <- terra::rast(system.file("ex/elev.tif", package="terra"))

# 原始范围与扩展范围
ext_original <- terra::ext(r)
ext_expanded <- terra::ext(
  xmin(ext_original) - 1,
  xmax(ext_original) + 1,
  ymin(ext_original) - 1,
  ymax(ext_original) + 1
)

# 下载对应扩展范围的OSM瓦片
osm_rast <- maptiles::get_tiles(ext_expanded, provider = "OpenStreetMap", crop = TRUE, zoom = 8)

# 重采样栅格以匹配瓦片分辨率
transcura_resamp <- resample(r, osm_rast, method = 'bilinear')
plot(osm_rast, axes = TRUE, mar = c(3.1, 3.1, 2.1, 7.1))
plot(transcura_resamp, add = TRUE, axes = TRUE, grid = TRUE)

但当调用add_grid()函数添加自定义网格后，网格出现错位异常：

library(terra)
library(maptiles)

# 示例栅格：局部高程数据
r <- terra::rast(system.file("ex/elev.tif", package="terra"))

# 原始范围与扩展范围
ext_original <- terra::ext(r)
ext_expanded <- terra::ext(
  xmin(ext_original) - 1,
  xmax(ext_original) + 1,
  ymin(ext_original) - 1,
  ymax(ext_original) + 1
)

# 下载对应扩展范围的OSM瓦片
osm_rast <- maptiles::get_tiles(ext_expanded, provider = "OpenStreetMap", crop = TRUE, zoom = 8)

# 重采样栅格以匹配瓦片分辨率
transcura_resamp <- resample(r, osm_rast, method = 'bilinear')
plot(osm_rast, axes = TRUE, mar = c(3.1, 3.1, 2.1, 7.1))
plot(transcura_resamp, add = TRUE, axes = TRUE, grid = TRUE)
add_grid(col = "black", lty = "dotted")

原因分析 #

add_grid()默认会基于当前绘图设备的默认坐标系统生成网格，但OSM瓦片采用的是Web Mercator投影（EPSG:3857），而原始栅格r是WGS84地理坐标系（EPSG:4326）。虽然通过resample()让栅格匹配了瓦片的分辨率，但add_grid()没有自动识别当前绘图的投影范围，导致网格错位。

解决方法 #

方法1：直接指定栅格对象给add_grid() #

调用add_grid()时传入当前绘图的栅格（如osm_rast），让函数基于该栅格的范围和投影生成网格：

# 替换原add_grid行
add_grid(osm_rast, col = "black", lty = "dotted")

方法2：提前统一投影 #

先将原始栅格转换为Web Mercator投影，再进行后续操作，确保全流程坐标系统一致：

library(terra)
library(maptiles)

r <- terra::rast(system.file("ex/elev.tif", package="terra"))
# 转换原始栅格到Web Mercator投影（与OSM瓦片一致）
r_3857 <- project(r, "EPSG:3857")

ext_original <- terra::ext(r_3857)
ext_expanded <- terra::ext(
  xmin(ext_original) - 1000,  # 投影单位为米，调整合适的扩展距离
  xmax(ext_original) + 1000,
  ymin(ext_original) - 1000,
  ymax(ext_original) + 1000
)

osm_rast <- maptiles::get_tiles(ext_expanded, provider = "OpenStreetMap", crop = TRUE, zoom = 8)
transcura_resamp <- resample(r_3857, osm_rast, method = 'bilinear')

plot(osm_rast, axes = TRUE, mar = c(3.1, 3.1, 2.1, 7.1))
plot(transcura_resamp, add = TRUE, axes = TRUE, grid = TRUE)
add_grid(col = "black", lty = "dotted")

两种方法都能让生成的网格与初始正常网格一致。

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