Блог

У меня есть тысячи частично перекрывающихся растров, каждый размером около 6500 x 6500 пикселей, с перекрытием 500 пикселей на каждом из четырех краев.

Я пытаюсь объединить их в один геотиф и усреднить перекрывающиеся области. Я следовал указаниям в этом вопросе, которые привели к следующим выводам.

1. gdal_merge и gdal_translate -r average варианты на самом деле не усредняют перекрывающиеся регионы
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))
 

Это работает, но это:

а) Загружает все TIF в оперативную память-это около 30 ГБ, что грубо б) Для запуска требуется несколько дней

Любопытно, знает ли кто-нибудь более быстрый способ решения этой проблемы!