使用gdaltranslate.exe处理重叠栅格的困惑与最大值提取需求

为什么你的NA值没保留？ #

这其实是gdalbuildvrt的默认行为导致的：它会按照输入文件的顺序，优先取第一个非NA的像素值。也就是说，当多个栅格叠加时，只要前面的栅格在某个位置有有效值，就会直接用这个值，完全忽略后面栅格在该位置的内容（包括NA）。

看你的第一个例子：

• file_a.tif的最后5个像素是96-100（全是有效值）
• file_b.tif的最后5个像素是NA
  所以VRT会直接用file_a的有效值，自然不会保留file_b的NA。

至于你调换顺序后得到全1-100的结果，应该是你把file_b.tif放在了输入列表的第一位？因为如果file_b在前，它的全有效值会覆盖后面file_a的101-200，所以最终结果就是file_b的值。

如何实现提取重叠栅格的最大值？ #

既然你的最终目标是对多个栅格取像素最大值，默认的GDAL叠加逻辑肯定不满足需求，这里给你两个最实用的方案：

方案1：用terra直接处理（最适合你的R环境）

你已经在使用R的terra包，完全可以跳过GDAL命令行，直接在R里完成最大值计算，既直观又不容易出错：

library(terra)

# 读取两个栅格文件
a <- rast("file_a.tif")
b <- rast("file_b.tif")

# 计算像素级最大值，自动处理NA（na.rm=TRUE表示忽略NA取有效最大值）
max_raster <- max(a, b, na.rm = TRUE)

# 写入结果，带上你需要的压缩参数
writeRaster(max_raster, "max_result.tif",
            options = c("PREDICTOR=2", "COMPRESS=LZW", "TILED=YES"))

# 验证结果
minmax(max_raster)

这个方法会对每个像素取两个栅格中的最大值——如果其中一个是NA，就取另一个的有效值；如果两个都是有效值，就取更大的那个，完全符合你的需求。

方案2：用GDAL命令行实现最大值计算

如果你一定要用GDAL工具链，可以用gdal_calc.py来完成像素级的最大值计算，步骤如下：

1. 先把多个栅格合并成一个多波段VRT（每个输入栅格对应一个波段）：

"C://Program Files//QGIS 3.26.1//bin//gdalbuildvrt.exe" -separate -input_file_list vrt_list.txt temp_multi_band.vrt

2. 用gdal_calc.py对多波段VRT的所有波段取像素最大值：

"C://Program Files//QGIS 3.26.1//bin//python.exe" "C://Program Files//QGIS 3.26.1//bin//gdal_calc.py" -A temp_multi_band.vrt --outfile max_result.tif --calc "numpy.max(A, axis=0)" --co "COMPRESS=LZW" --co "TILED=YES" --co "PREDICTOR=2"

这个命令会调用numpy对每个像素位置的所有波段值取最大值，自动处理NA值，输出符合要求的TIFF文件。

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