平均数千个与GDAL部分重叠的光栅

Question

我有数千个部分重叠的光栅，每个大约6500 x 6500 px，四条边上各有500 px的重叠

我试图将它们合并到一个geotiff中，并平均重叠区域。我遵循了问题this中的指导，得出了以下结论

1. gdal_merge和gdal_translate{}选项实际上并不平均重叠区域

2. 使用像素函数和GDALBuildVRT可以工作

基于此，我设置了以下代码来创建VRT、添加像素函数和创建单个TIF

out_name = "output"
tifs = ['tif1.tif', .... 'tif1235.tif']

gdal.BuildVRT(f'{out_name}.vrt', tifs, options=gdal.BuildVRTOptions(srcNodata=255, VRTNodata=255))

add_pixel_fn(f'{out_name}.vrt')

ds = gdal.Open(f'{out_name}.vrt')
translateoptions = gdal.TranslateOptions(gdal.ParseCommandLine("-ot Byte -co COMPRESS=LZW -a_nodata 255"))
ds = gdal.Translate(f'{out_name}.tif', ds, options=translateoptions)

其中，像素函数定义为

def add_pixel_fn(filename: str) -> None:
    """inserts pixel-function into vrt file named 'filename'
    Args:
        filename (:obj:`string`): name of file, into which the function will be inserted
        resample_name (:obj:`string`): name of resampling method
    """

    header = """  <VRTRasterBand dataType="Byte" band="1" subClass="VRTDerivedRasterBand">"""
    contents = """
    <PixelFunctionType>average</PixelFunctionType>
    <PixelFunctionLanguage>Python</PixelFunctionLanguage>
    <PixelFunctionCode><![CDATA[
from numba import jit
import numpy as np
@jit(nogil=True)
def average_jit(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize, raster_ysize, buf_radius, gt):
    np.mean(in_ar, axis = 0,out = out_ar, dtype = 'uint8')

def average(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,raster_ysize, buf_radius, gt):
    average_jit(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,raster_ysize, buf_radius, gt)
]]>
    </PixelFunctionCode>"""

    lines = open(filename, 'r').readlines()
    lines[3] = header
    lines.insert(4, contents)
    open(filename, 'w').write("".join(lines))

这是可行的，但它：

a）将所有TIF加载到RAM中——这大约是30+gb，这很粗糙 b） 运行需要几天时间

好奇是否有人知道任何更快的方法来解决这个问题