#python #raster #gdal #rasterio
Вопрос:
У меня есть тысячи частично перекрывающихся растров, каждый размером около 6500 x 6500 пикселей, с перекрытием 500 пикселей на каждом из четырех краев.
Я пытаюсь объединить их в один геотиф и усреднить перекрывающиеся области. Я следовал указаниям в этом вопросе, которые привели к следующим выводам.
gdal_merge
иgdal_translate
-r average
варианты на самом деле не усредняют перекрывающиеся регионы- Использование функции пикселей и 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 ГБ, что грубо б) Для запуска требуется несколько дней
Любопытно, знает ли кто-нибудь более быстрый способ решения этой проблемы!
Комментарии:
1. Я удивлен, что он пытается загрузить все файлы в память, попробуйте написать плиточный TIFF (
-co TILED=YES
)