Объединение перекрывающихся растров по минимальному значению в python

Question

Я хотел бы объединить некоторые 1-полосные растровые изображения. Я бы хотел, чтобы перекрытие обрабатывалось так, чтобы было выбрано наименьшее значение в пикселе. Все изображения имеют одинаковую проекцию.

Я попытался посмотреть на gdalwarp и gdal_merge (в командной строке), но в перекрытии они просто используют значения из последнего изображения.

Я также видел предложения использовать смесь PIL или пасту, но для этого требуется альфа-слой, который указывает, как следует смешивать перекрытия, чего, конечно, у меня нет.

Есть ли кто-нибудь, кто знает, как сделать перекрытие зависимым от фактических значений на изображениях?

Answer 1

Насколько я знаю, нет способа сделать это с помощью стандартных утилит командной строки. Вы попробуйте 'gdal_calc.py', который обычно поставляется с GDAL, для этого требуется, чтобы растр имел одинаковые размеры, поэтому, возможно, потребуется некоторая предварительная обработка (например, с gdalwarp). Минимум двух растров можно рассчитать с помощью чего-то вроде:

python gdal_calc.py -A file1.tif -B file2.tif --calc="minimum(A,B)" --outfile=res.tif

Я не совсем уверен, что происходит в регионах без наложения, возможно, вам нужно добавить ключевое слово nodata, более подробную информацию можно получить по адресу: http://www.gdal.org/gdal_calc.html

Вы должны указать полный путь к gdal_calc.py,

Для будущих выпусков:

Начиная с GDAL 2.2, есть поддержка функций VRT Pixel, написанных на Python:

http://www.gdal.org/gdal_vrttut.html

После этого вы сможете сделать VRT похожим на:

<VRTDataset rasterXSize="20" rasterYSize="20">
  <SRS>EPSG:26711</SRS>
  <GeoTransform>440720,60,0,3751320,0,-60</GeoTransform>
  <VRTRasterBand dataType="Byte" band="1" subClass="VRTDerivedRasterBand">
    <PixelFunctionType>add</PixelFunctionType>
    <PixelFunctionLanguage>Python</PixelFunctionLanguage>
    <PixelFunctionCode><![CDATA[
import numpy as np
def add(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,
                   raster_ysize, buf_radius, gt, **kwargs):
    np.round_(np.clip(np.minimum(in_ar, axis = 0, dtype = 'uint16'),0,255),
              out = out_ar)
]]>
    </PixelFunctionCode>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">byte.tif</SourceFilename>
    </SimpleSource>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">byte2.tif</SourceFilename>
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

Кроме того, начиная с GDAL 2.2, будут встроены функции пикселей, что является отличной идеей. Я пока не вижу агрегации, такие как min, max, означает, но это, вероятно, будет легко добавить, когда инфраструктура будет там. http://www.gdal.org/gdal_vrttut.html