python-2.7 - 使用来自另一个栅格的数据填充 geotiff 中的无数据值

Question

我有两个不同分辨率的重叠光栅文件(A=25m 和 B=5m)。我想以最快最有效的方式(如果可能使用 python 绑定(bind)进行自动化)用 A 中的值填充 B 中的无数据值 (-999)。

我测试了不同的方法:

1. 使用 gdalwarp 将光栅 A 重采样为 B 的分辨率，然后使用 gdal_calc 计算公式 gdal_calc.py -A MNT_5m.tif -B BATI_5m.tif --outfile=out.tif --calc="maximum( A,B)”，但这并没有做任何事情，而且重新采样是一项耗时的操作。
2. 再次使用 gdal_calc 计算公式 --calc="A*(B=-999)+B*(B!=-999)" 但这给了我以下错误:SyntaxError:语法无效。

谁能帮我解决这个问题？

Answer 1

如果 nodata 值设置正确(在元数据中)，使用 gdal.BuildVRT 可能是执行此操作的好方法。

BuildVRT 允许分辨率差异，因此您或许可以跳过 warp 步骤。但是所有文件的投影必须相同。

命令中首先列出的文件最终“位于顶部”，因此应该是您的 5m 文件，然后是 25m 文件。

import gdal

dsvrt = gdal.BuildVRT('output.vrt', ['MNT_5m.tif', 'BATI_25m.tif'],
                      resolution='highest', resampleAlg=gdal.GRA_Bilinear)

如果作为输出的 VRT 文件不够，您可以轻松地使用 Translate 将其转换为 TIFF 或其他格式:

dstif = gdal.Translate('output.tif', dsvrt)

dsvrt = None
dstif = None

如果确实需要使用gdal.Warp，可以在gdal.BuildVRT之前添加，并将返回的Dataset传入BuildVRT (与作为字符串的文件名相反)。

这种工作方式的一个好处是您可以轻松地在是否输出中间文件之间切换。通过使用 /vsimem/tempfile 作为输出，并将 VRT 作为输出格式，您可以即时(在内存中)执行所有中间步骤。如果您想要输出，用于调试或其他目的，将文件名更改为磁盘上的某个位置即可获得中间结果。