Подводные камни gdal_calc

Не так давно я писал о растровом калькуляторе в GDAL. Вещь это нужная и обладающая достаточным функционалом для решения большинства задач. Но, как оказалось, не лишенная багов.

Баловался с этим калькулятором на досуге, заглянул в код. Насторожило то, что данные читаются «как есть», т.к. когда-то я на эти грабли уже наступал (/me на несколько минут погружается в воспоминания. Веселое было время: два параллельных проекта, оба связаны с GDAL; куча нового материала…). Небольшая проверка подтвердила мои опасения. Вот результат расчета NDVI для растра с типом данных Byte в RasterCalc (сохранен как Float32)

А вот это результат расчета при помощи gdal_calc.py

Картинкой это не передашь, но белые области на самом деле не белые, а содержат значение NODATA.

Что и говорить, результаты очень интересные. Но объясняются они просто: т.к. чтение данных происходит «как есть», то при выполнении операций  можно нарваться на такие вещи как целочисленное переполнение и округление результата (в Python результат деления двух целых тоже целое число).

Ошибка налицо. Пишем автору, что он козелсообщаем о баге и предлагаем патч.

Использующие этот инструмент могут не ждать выхода следующей версии GDAL и исправить ошибку самостоятельно. Для этого открываем gdal_calc.py в любимом текстовом редакторе, и строки 230-232

myval=BandReadAsArray(myFiles[i].GetRasterBand(myBands[i]), xoff=myX, yoff=myY, win_xsize=nXValid, win_ysize=nYValid)

приводим к виду

myval=BandReadAsArray(myFiles[i].GetRasterBand(myBands[i]), xoff=myX, yoff=myY, win_xsize=nXValid, win_ysize=nYValid).astype(numpy.float32)
Мітки:

Залишити відповідь

Ваша e-mail адреса не оприлюднюватиметься. Обов’язкові поля позначені *

*