【GDAL】如何利用GDAL批量进行等同于Arcmap中Mosaic To New Raster的操作

GDAL.BuildVRT

这里主要参考 gdal官方文件，生成VRT格式文件：GDAL.BuildVRT

# 最好保证所有栅格数据在同一坐标系下，且具有相同的分辨率
vrt = gdal.BuildVRT(f"{name}.vrt", tif_files_list) # tif_files_list 包含所有要Mosaic的栅格文件
vrt = None

VRT Pixelfunction

随后利用 lxml库读写VRT文件，并加入自定义的Pixel function以解决重叠栅格取值的问题，下面代码以获取最大值为例。如果对VRT文件的构成感兴趣可以参考: VRT

from lxml import etree


vrt_tree = etree.parse(f"{name}.vrt")
vrt_root = vrt_tree.getroot()
vrtband1 = vrt_root.findall(".//VRTRasterBand[@band='1']")[0]

vrtband1.set("subClass", "VRTDerivedRasterBand")
pixelFunctionType = etree.SubElement(vrtband1, 'PixelFunctionType')
pixelFunctionType.text = "find_max"
pixelFunctionLanguage = etree.SubElement(vrtband1, 'PixelFunctionLanguage')
pixelFunctionLanguage.text = "Python"
pixelFunctionCode = etree.SubElement(vrtband1, 'PixelFunctionCode')
pixelFunctionCode.text = etree.CDATA("""
import numpy as np

def find_max(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize, raster_ysize, buf_radius, gt, **kwargs):
	np.amin(in_ar, axis=0, initial=255, out=out_ar) 
""")

GDAL.Translate

最后利用gdal.Translate将VRT文件转换为TIF格式，同时可以利用该函数对TIF文件进行压缩：gdal.Translate

gdal.Translate(f"{name}.TIF", f"{name}.vrt")