Растровая арифметика в GDAL

Растровый калькулятор это круто, но когда возникает необходимость обработать кучу растров начинаются проблемы. Ну нет в калькуляторе (ни в одном, ни в другом) пакетного режима, нет. Вот в моем TODO он есть. Правда, от осознания этого факта легче не становится, когда у тебя 100500 растров для которых надо расcчитать, к примеру NDVI или выполнить замену пикселей по хитрому условию.

Тут на помощь приходит GDAL и Python. Буквально за 10 минут пишется скрипт, благо есть документация и на английском, и на русском. Но опять же, писать скрипты на каждый чих как-то не очень целесообразно.

Открою маленький секрет. В GDAL 1.8.0 есть замечательный инструмент с незатейливым названием gdal_calc.py. Этот небольшой (около 300 строк) скрипт на Python работает с растрами, имеющими одинаковые размеры (проверка на соответствие проекций не выполняется), и поддерживает базовые арифметические и логические действия. Пользоваться просто:

# сумма двух растров
gdal_calc.py -A input1.tif -B input2.tif --outfile=result.tif --calc="A+B"
# среднее значение двух растров
gdal_calc.py -A input1.tif -B input2.tif --outfile=result.tif --calc="(A+B)/2"
# разность каналов
gdal_calc.py -A input.tif --A_band=1 -B input.tif --B_band=2 --outfile=result.tif --calc="A-B"

На вход можно подавать до 26 растров, этого должно хватить для большинства применений. Как видно из примеров, есть поддержка скобок, можно манипулировать каналами, а консольная природа позволяет легко использовать скрипт для пакетной обработки. Его бы чуть-чуть допилить и получится отличный консольный калькулятор.