基于栅格单元键值实现表格与栅格的左连接（大栅格适配）

需求：以栅格的单元值为键，将表格数据关联到栅格，需适配1亿至100亿单元的超大规模栅格，当前方案存在内存瓶颈，寻求terra包内置功能或data.table优化技巧。

当前解决方案（存在内存限制） #

# Load Libs
library(data.table)
library(terra)

# -------- example data ------------

# Create Raster 
r <- rast(matrix(c(1:24), nrow = 4, ncol = 6))
names(r) <- "r_key"

# Create tabular data
dt <- data.table(
  r_key = c(1:50), 
  vals = sample(c(1000:5000), 50)
)

# -------- end example data ----------

# -------- Current Solution ----------

# 该方案在栅格单元数增长时会遇到内存瓶颈
 
r_table <- data.table()
r_table[, r_key := values(r)]
# 创建索引以保证栅格单元格顺序
r_table[, index := .I]

# 常规合并操作
r_table <- merge(r_table,
                 dt,
                 by = "r_key",
                 all.x = TRUE)

# 恢复原始顺序
setorder(r_table, "index")

# 生成带新值的栅格
r_new <- deepcopy(r)
values(r_new) <- r_table[["vals"]]

plot(r_new)

# ---------- end current solution ------

高效优化方案 #

1. 使用terra包内置的classify函数（推荐） #

terra的classify函数专门用于栅格值的映射替换，底层基于向量化操作，无需创建全量单元格数据表，内存占用极低，完美适配超大规模栅格场景。

只需将表格数据转换为两列规则矩阵（第一列为栅格键值，第二列为对应替换值），直接传入函数即可：

# 转换规则格式：r_key -> vals
rules <- as.matrix(dt)
# 执行键值映射，na.rm=FALSE保留未匹配的NA值
r_new <- classify(r, rules, na.rm = FALSE)
plot(r_new)

2. data.table向量级查询优化 #

跳过中间数据表的创建，直接对栅格值向量做键值匹配，避免merge和排序步骤，大幅降低内存消耗：

# 提取栅格所有单元值作为键向量
r_keys <- values(r)
# 用data.table快速匹配，返回对应vals，未匹配项自动为NA
new_vals <- dt[J(r_keys), vals, on = "r_key"]
# 赋值给新栅格
r_new <- deepcopy(r)
values(r_new) <- new_vals
plot(r_new)

这种方式仅处理两个向量（r_keys和new_vals），无需存储带索引的全量数据表，内存效率远高于原方案。

3. 百亿级栅格的分块处理 #

如果栅格规模大到无法一次性加载到内存，使用terra的分块读写API，逐块完成匹配写入：

# 创建与原栅格结构一致的空输出栅格
r_new <- rast(r)
# 初始化分块写入（指定输出文件路径）
writeStart(r_new, filename = "output_raster.tif", overwrite = TRUE)
# 按每块100万单元设置分块数量（可根据内存调整）
block_size <- 1000000
n_blocks <- ceiling(ncell(r) / block_size)

for (i in 1:n_blocks) {
  # 计算当前块的起始位置和单元数
  start_pos <- (i-1)*block_size + 1
  n_cells <- min(block_size, ncell(r) - start_pos + 1)
  # 读取当前块的单元值
  block_vals <- readValues(r, start = start_pos, n = n_cells)
  # 匹配表格数据
  matched_vals <- dt[J(block_vals), vals, on = "r_key"]
  # 写入当前块
  writeValues(r_new, matched_vals, start = start_pos)
}
# 结束分块写入
writeStop(r_new)

通过分块读取-匹配-写入的方式，避免一次性加载百亿级单元数据，彻底解决内存溢出问题。

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