elevatr全球高程数据转DataFrame遇内存不足，求助terra blocks用法

问题场景 #

我用elevatr包获取全球高程数据，耗时约2小时，保存后用terra::rast按需加载：

world_elevation <- elevatr::get_elev_raster(
     locations = world_poly,
     z = 7,
     clip = "locations",
     override_size_check = T
)

writeRaster(world_elevation, 
            "E:\\Research\\Study\\calvarium\\Maps\\globElev\\world_elevation.tif", 
            options = c("COMPRESS = NONE", "TFW = YES"))
globRaster <- terra::rast("E:\\Research\\Study\\calvarium\\Maps\\globElev\\world_elevation.tif")

尝试将栅格转为带坐标的数据框时，因16G内存不足报错：

globRaster_df <- globRaster |>
    as.data.frame(xy = T) |>
    na.omit()

错误信息：

Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'object' in selecting a method for function 'na.omit': cannot allocate vector of size 23.6 Gb.

我希望获取完整的globRaster_df，注意到terra的blocks函数可分块读取栅格，但不清楚具体用法。后续计划用以下代码绘图：

globRaster_df <- country_elevation |>
    as.data.frame(xy = T) |>
    na.omit()

names(country_elevation_df)[3] <- "elevation"

country_map <-
    ggplot(data = globRaster_df ) +
    geom_raster(
        aes(x = x, y = y, fill = elevation),
        alpha = 1
    ) +
    coord_sf() +
    labs(
        x = "",
        y = "",
        title = "",
        subtitle = "",
        caption = ""
    ) +
    theme_minimal() +
    theme(
        axis.line = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.margin = unit(c(t = 0, r = 0, b = 0, l = 0), "cm"),
        plot.background = element_blank(),
        panel.background = element_blank(),
        panel.border = element_blank()
    )

解决方案 #

方案1：分块读取栅格并合并数据框 #

利用terra::blocks获取栅格的分块信息，循环读取每个块并转换为数据框，最后合并，避免一次性加载全部数据：

library(terra)
library(dplyr)
library(tidyr)

# 获取栅格的分块参数
b <- blocks(globRaster)
df_list <- list()

# 遍历每个块处理
for (i in 1:nrow(b)) {
  # 启动栅格读取会话
  rast_obj <- readStart(globRaster)
  # 读取当前块的数值
  block_vals <- readValues(rast_obj, 
                           row = b$row[i], nrows = b$nrows[i],
                           col = b$col[i], ncols = b$ncols[i])
  readStop(rast_obj)
  
  # 生成当前块的坐标网格并绑定高程值
  block_df <- expand_grid(
    x = seq(b$xmin[i], b$xmax[i], length.out = b$ncols[i]),
    y = seq(b$ymax[i], b$ymin[i], length.out = b$nrows[i])
  ) |>
    mutate(elevation = as.vector(block_vals)) |>
    na.omit()
  
  df_list[[i]] <- block_df
}

# 合并所有块的数据框
globRaster_df <- bind_rows(df_list)

方案2：直接用SpatRaster绘图（更高效） #

无需转换为超大数据框，ggplot2配合terra可以直接读取磁盘上的栅格文件分块绘图，大幅节省内存：

library(ggplot2)
library(terra)

country_map <- ggplot() +
  # 直接传入SpatRaster对象
  geom_spatraster(data = globRaster, aes(fill = elevation)) +
  coord_sf() +
  labs(x = "", y = "", title = "", subtitle = "", caption = "") +
  theme_minimal() +
  theme(
    axis.line = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    legend.position = "none",
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    plot.margin = unit(c(t = 0, r = 0, b = 0, l = 0), "cm"),
    plot.background = element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank()
  )

# 保存图像
ggsave("global_elevation_map.png", country_map, width = 12, height = 6, dpi = 300)

推荐方案2，因为它不需要将整个栅格数据加载到内存，直接从磁盘分块读取绘图，适合处理大尺寸栅格。

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