如何在terra::mosaic()中使用自定义众数函数？大GeoTIFF拼接方案咨询

一、R terra包自定义众数函数实现 #

terra的mosaic()函数支持传入自定义fun参数，你可以自己编写众数计算函数来处理重叠像素。针对大数据量，优先用虚拟栅格（VRT）减少内存占用。

自定义众数函数与拼接代码 #

# 高效众数计算函数（自动过滤NA值）
mode_fun <- function(x, na.rm = TRUE) {
  if (na.rm) x <- x[!is.na(x)]
  if (length(x) == 0) return(NA)
  tab <- table(x)
  as.numeric(names(tab)[which.max(tab)])
}

library(terra)
# 读取所有影像路径
files <- list.files(path = "你的影像目录", pattern = "\\.tif$", full.names = TRUE)
# 创建虚拟栅格（无需一次性加载所有影像到内存）
vrt <- vrt(files)
# 执行众数拼接
mosaic_result <- mosaic(vrt, fun = mode_fun)
# 保存结果
writeRaster(mosaic_result, "mosaic_mode.tif", overwrite = TRUE)

大数据量优化 #

• 用vrt()替代直接读取所有栅格，terra会按需加载数据，大幅降低内存压力。
• 可通过terraOptions(memfrac = 0.5)限制terra使用的内存比例，避免内存溢出。

二、GDAL工具替代方案 #

针对6000张影像的规模，GDAL命令行+自定义Python脚本的组合性能更优，推荐以下两种方式：

方式1：rasterio分块处理 #

先创建虚拟栅格，再分块读取计算众数：

# 第一步：创建所有影像的虚拟栅格
gdalbuildvrt input_images.vrt /path/to/your/images/*.tif

import rasterio
import numpy as np
from rasterio.merge import merge

def calculate_mode(arr):
    mask = np.isnan(arr) if arr.dtype in [np.float32, np.float64] else (arr == -9999)
    arr = arr[~mask]
    if len(arr) == 0:
        return np.nan
    vals, counts = np.unique(arr, return_counts=True)
    return vals[np.argmax(counts)]

# 读取虚拟栅格
src = rasterio.open('input_images.vrt')
# 定义输出参数（可根据需求调整尺寸和分辨率）
out_width, out_height = 10000, 10000
out_transform = src.transform * src.transform.scale(
    src.width / out_width,
    src.height / out_height
)

# 分块写入结果
with rasterio.open(
    'mosaic_mode_gdal.tif',
    'w',
    driver='GTiff',
    height=out_height,
    width=out_width,
    count=src.count,
    dtype=src.dtypes[0],
    crs=src.crs,
    transform=out_transform,
    nodata=src.nodata
) as dst:
    for _, window in dst.block_windows(1):
        # 读取当前块的所有影像数据
        arr = src.read(window=window, boundless=True)
        # 对每个像素维度计算众数
        mode_arr = np.apply_along_axis(calculate_mode, 0, arr)
        dst.write(mode_arr, window=window)

方式2：gdal_grid（适合离散型数据） #

将影像像素转成点要素后，用gdal_grid的众数插值生成结果：

# 循环将所有影像转成XYZ格式点文件（示例）
for file in /path/to/your/images/*.tif; do
  gdal_translate -of XYZ "$file" "${file%.tif}.txt"
done
# 合并所有点文件
cat /path/to/your/images/*.txt > all_points.txt
# 生成众数拼接影像
gdal_grid -a mode:radius=1 -txe xmin xmax -tye ymax ymin -outsize 10000 10000 -of GTiff all_points.txt mosaic_mode_grid.tif

三、方案对比 #

• R terra：代码简洁，适合熟悉R环境的用户，VRT模式可应对大部分大规模场景。
• GDAL/Python：底层读写效率更高，分块处理机制更适合超大规模数据，性能优势明显。

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