Растровая арифметика в QGIS с помощью RasterCalc

[updates-bot]

Обсуждение статьи "Растровая арифметика в QGIS с помощью RasterCalc"

http://gis-lab.info/qa/rastercalc.html

[kryl]

Не хочет работать с растром после атмосферной коррекции (значения каналов Float)

Код: Выделить всё

Ошибка при выполнении сценария Python:

Traceback (most recent call last):
  File "C:/Users/ws_3/.qgis/python/plugins\rastercalc\rastercalc.py", line 160, in run
    dlg = rastercalcdialog.RasterCalcDialog()
  File "C:/Users/ws_3/.qgis/python/plugins\rastercalc\rastercalcdialog.py", line 106, in __init__
    self.manageGui()
  File "C:/Users/ws_3/.qgis/python/plugins\rastercalc\rastercalcdialog.py", line 114, in manageGui
    self.layerNames = [ str( z.name() ) for z in self.mapLayers ]
UnicodeEncodeError: 'ascii' codec can't encode character u'\u0441' in position 21: ordinal not in range(128)

Версия Python:
2.5.1 (r251:54863, Apr 18 2007, 08:51:08) [MSC v.1310 32 bit (Intel)]


Версия QGIS:
1.6.0-Trunk Trunk, 14423

Путь поиска Python: ['C:/OSGeo4W/apps/qgis-dev/./python', 'C:/Users/ws_3/.qgis/python', 'C:/Users/ws_3/.qgis/python/plugins', 'C:/OSGeo4W/apps/qgis-dev/./python/plugins', 'C:\\OSGeo4W\\apps\\gdal-17\\pymod', 'C:\\', 'C:\\Windows\\system32\\python25.zip', 'C:\\OSGeo4W\\apps\\Python25\\DLLs', 'C:\\OSGeo4W\\apps\\Python25\\lib', 'C:\\OSGeo4W\\apps\\Python25\\lib\\plat-win', 'C:\\OSGeo4W\\apps\\Python25\\lib\\lib-tk', 'C:\\OSGeo4W\\apps\\qgis-dev\\bin', 'C:\\OSGeo4W\\apps\\Python25', 'C:\\OSGeo4W\\apps\\Python25\\lib\\site-packages', 'C:\\OSGeo4W\\apps\\Python25\\lib\\site-packages\\win32', 'C:\\OSGeo4W\\apps\\Python25\\lib\\site-packages\\win32\\lib', 'C:\\OSGeo4W\\apps\\Python25\\lib\\site-packages\\Pythonwin', 'C:\\OSGeo4W\\apps\\Python25\\lib\\site-packages\\wx-2.8-msw-unicode', 'C:\\Users\\ws_3\\.qgis\\python\\plugins\\cadtools\\tools', 'C:\\OSGeo4W\\apps\\qgis-dev\\python\\plugins\\fTools\\tools', '/usr/share/qgis/python', '/usr/share/qgis/python', '.', 'C:\\Users\\ws_3\\.qgis\\python\\plugins\\ziplayers\\logic', 'C:\\Users\\ws_3\\.qgis\\python\\plugins\\ziplayers\\gui', '.']

Инфа о растре (слегка сокращенная):

Код: Выделить всё

Driver: HFA/Erdas Imagine Images (.img)
Files: D:/arcwork/yaroslavl2/space/l5cor/LT51780192007229MOR00/LT51780192007229MOR00ñ.img
Size is 8581, 7732
Coordinate System is:
PROJCS["UTM Zone 37, Northern Hemisphere",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",39],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1],
    AUTHORITY["EPSG","32637"]]
Origin = (413679.767395963720000,6621345.494110753800000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Corner Coordinates:
Upper Left  (  413679.767, 6621345.494) ( 37d27'54.84"E, 59d43'15.82"N)
Lower Left  (  413679.767, 6389385.494) ( 37d33'14.85"E, 57d38'18.63"N)
Upper Right (  671109.767, 6621345.494) ( 42d 2'25.65"E, 59d41'41.34"N)
Lower Right (  671109.767, 6389385.494) ( 41d51'52.51"E, 57d36'51.54"N)
Center      (  542394.767, 6505365.494) ( 39d43'52.62"E, 58d41'11.07"N)
Band 1 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = Layer_1
  Min=0.000 Max=0.342 
  Minimum=0.000, Maximum=0.342, Mean=0.098, StdDev=0.010
  Overviews: 4291x3866, 2146x1933, 1073x967, 537x484, 269x242, 135x121, 68x61
  Metadata:
    STATISTICS_MINIMUM=3.3208416425623e-006
    STATISTICS_MAXIMUM=0.34183913469315
    STATISTICS_MEAN=0.097974001067624
    STATISTICS_MEDIAN=0.09480934733952
    STATISTICS_MODE=0.092138755043806
    STATISTICS_STDDEV=0.010419113221324
    STATISTICS_HISTONUMBINS=256
    STATISTICS_HISTOMIN=3.3208416425623e-006
    STATISTICS_HISTOMAX=0.34183913469315
    STATISTICS_SKIPFACTORX=2
    STATISTICS_SKIPFACTORY=2
    STATISTICS_EXCLUDEDVALUES=0
    LAYER_TYPE=athematic
    STATISTICS_HISTOBINVALUES=10736|1504|1040|916|644|680|572|536|564|628|708|560|588|564|644|780|2544|14160|17764|9616|5332|3424|2904|2196|1908|1736|1452|960|544|388|292|716|2012|2740|6900|11272|18020|16264|10960|9204|6404|4020|2192|2376|2232|2300|2200|2532|2996|2772|2536|3328|5252|6624|7392|7624|6192|4928|4052|3456|2576|2560|15528|123368|474556|1050524|1620932|2156696|2801436|3090948|2962460|2691028|2445744|2167564|1752536|1339576|1084976|947172|868868|817836|729064|582868|464828|382572|221532|288648|246340|216352|193832|175524|160572|142724|125320|111028|101336|93340|81320|67728|56876|50204|46100|41916|34268|26904|21316|18404|16556|14200|10924|7972|5960|4896|4072|3464|2364|1608|1228|980|728|588|452|384|340|244|260|240|232|100|132|76|76|52|56|32|24|16|32|4|16|4|4|8|8|12|8|0|8|12|4|8|4|8|0|0|8|0|8|0|12|4|8|0|0|4|0|0|0|4|0|8|0|4|4|0|0|0|0|0|0|0|12|0|0|0|0|0|0|4|0|0|8|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|4|0|0|0|0|0|0|0|4|0|0|4|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|4|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|4|
  Image Structure Metadata:
    COMPRESSION=RLE


Band 2 Block=512x512 Type=Float32, ColorInterp=Undefined
  Description = Layer_2
  Min=0.000 Max=0.528 
  Minimum=0.000, Maximum=0.528, Mean=0.081, StdDev=0.013
  Overviews: 4291x3866, 2146x1933, 1073x967, 537x484, 269x242, 135x121, 68x61
  Metadata:
    STATISTICS_MINIMUM=2.9177457690821e-005
    STATISTICS_MAXIMUM=0.52833443880081
    STATISTICS_MEAN=0.081114446210786
    STATISTICS_MEDIAN=0.076385797261189
    STATISTICS_MODE=0.070194719979824
    STATISTICS_STDDEV=0.013150569806751
    STATISTICS_HISTONUMBINS=256
    STATISTICS_HISTOMIN=2.9177457690821e-005
    STATISTICS_HISTOMAX=0.52833443880081
    STATISTICS_SKIPFACTORX=2
    STATISTICS_SKIPFACTORY=2
    STATISTICS_EXCLUDEDVALUES=0
    LAYER_TYPE=athematic
    STATISTICS_HISTOBINVALUES=9228|1028|560|932|668|832|624|736|8060|28760|7808|6468|4560|3324|3896|1612|5432|10144|17164|25224|8992|7252|6860|7912|10276|6784|7556|10448|9312|12396|12004|58932|913448|2404960|6169956|3602532|3004264|1957324|1502892|1706636|1082964|1106924|1022656|1466964|990464|1074624|953956|882576|995628|453008|341572|241636|203068|262784|159796|111656|88508|67160|81868|47596|39268|33700|29128|33448|14560|11960|9768|12524|6684|5648|5000|4404|4836|2844|2368|2008|1804|2112|1216|940|864|788|1064|564|440|484|404|400|220|180|128|176|160|132|84|88|112|60|40|40|60|48|40|52|28|36|36|32|20|40|20|56|28|28|32|20|36|40|28|32|16|20|28|20|28|44|16|20|20|16|20|20|32|24|8|16|8|8|16|16|16|12|4|32|8|28|16|12|12|36|8|16|8|16|12|16|12|20|12|20|4|0|20|12|24|20|8|4|16|20|8|20|8|24|8|24|4|12|16|16|12|12|12|28|12|12|8|8|20|8|20|20|16|28|20|8|4|20|4|8|20|20|20|16|24|4|20|4|12|20|12|12|16|8|12|4|4|16|4|4|4|4|8|0|8|4|4|8|4|8|8|8|4|0|4|4|4|4|0|12|0|24|4|4|12|8|16|20|16|8|0|0|8|0|8|8|
  Image Structure Metadata:
    COMPRESSION=RLE
<GDALRasterAttributeTable Row0Min="2.917745769082103e-005" BinSize="0.002071785338600475">
 

[kryl]

Разобрался. Дело в русской "с" в названии файла

Код: Выделить всё

LT51780192007229MOR00ñ.img

[Voltron]

Доступен RasterCalc 0.2. В этой версии появились операторы сравнения и условные операторы.Также исправлены старые ошибки и добавлены новые

Более подробно узнать о новых функциях можно в обновленной статье.

[Denis Rykov]

Очень удобный инструмент, спасибо! А планируется ли возможность добавления пакетной обработки растров в RasterCalc?

[Voltron]

Пожалуйста . Вот уж не думал, что после появления в QGIS «родного» калькулятора, RasterCalc будет нужен кому-то, кроме автора.

Возможность пакетной обработки не планировалась, но если она нужна — запланируем. Как я понял, интересует возможность применения одного и того же выражения к нескольким выбранным растрам, правильно?

[Denis Rykov]

Так точно. А можно как-то воспользоваться плагином в обход QGIS? То есть, например, вызвать его из консоли и передать на вход нужное выражение.

А инструмент и вправду очень удобный и умеет кое-что, чего нет во встроенном калькуляторе.

[Максим Дубинин]

о, нас и на stack exchange, круто.

Я поддерживаю что нужен консольный растеркальк и напоминаю, что калькулятор растров был одним из гсоковских заданий для гдал, (nb: деньги нашим программирующим студентам не нужны, пнем их еще раз за лень). Но тут нужно крепко подумать, так как сейчас и в ядре кугиса есть растеркальк.

[Voltron]

В данный момент обходного варианта не существует, но есть несколько скриптов для GDAL разной степени продвинутости, реализующих растровую арифметику. Кугисовский калькулятор вынесен в API (пока только для C++), если его вынесут и в Python, то можно будет писать скрипты, но они будут привязаны к наличию в системе QGIS.

Если и делать консольный калькулятор, то с минимальными зависимостями (в идеале, только GDAL). Мне кажется, что отдельный калькулятор вещь нужная и должна существовать не зависимо от калькуляторов в QGIS.

[Максим Дубинин]

Мне кажется, что отдельный калькулятор вещь нужная и должна существовать не зависимо от калькуляторов в QGIS.

так точно

[Loup]

На всякий случай:

Код: Выделить всё

Ошибка при выполнении сценария Python:
...
UnicodeEncodeError: 'ascii' codec can't encode characters ... in position ...: ordinal not in range(128)

Ошибка появляется при использовании русских символов не только в имени файла, но и в пути к нему, например "D:\Дешифрирование\..."

[bolotoved]

Если нужно из многоканального растра удалить пиксели по маскирующему растру с определенным значением, то правильной ли будет следующая команда:

Код: Выделить всё

ne ( [mask]@1, 1, [raster])

?

При выполнении данной команды вылезает ошибка:

Код: Выделить всё

Ошибка при выполнении сценария Python:

Traceback (most recent call last):
  File "C:/Documents and Settings/admin/.qgis/python/plugins\rastercalc\rastercalcdialog.py", line 270, in accept
    result = rastercalcengine.evaluateStack( tmpStack, iter, sizeX, blk )
  File "C:/Documents and Settings/admin/.qgis/python/plugins\rastercalc\rastercalcengine.py", line 229, in evaluateStack
    return func[ op ]( inRaster, compare, replace )
  File "C:/Documents and Settings/admin/.qgis/python/plugins\rastercalc\rastercalcengine.py", line 95, in not_equal
    numpy.putmask( raster, tmp, replace )
ValueError: setting an array element with a sequence.

Версия Python:
2.5.2 (r252:60911, Feb 21 2008, 13:11:45) [MSC v.1310 32 bit (Intel)]


Версия:
1.8.0-Trunk Trunk, 67dfd2e

Путь поиска Python: ['C:/OSGeo4W/apps/qgis-dev/./python', 'C:/Documents and Settings/admin/.qgis/python', 'C:/Documents and Settings/admin/.qgis/python/plugins', 'C:/OSGeo4W/apps/qgis-dev/./python/plugins', 'C:\\OSGeo4W\\bin\\python25.zip', 'C:\\OSGeo4W\\apps\\Python25\\DLLs', 'C:\\OSGeo4W\\apps\\Python25\\lib', 'C:\\OSGeo4W\\apps\\Python25\\lib\\plat-win', 'C:\\OSGeo4W\\apps\\Python25\\lib\\lib-tk', 'C:\\OSGeo4W\\apps\\qgis-dev\\bin', 'C:\\OSGeo4W\\apps\\Python25', 'C:\\OSGeo4W\\apps\\Python25\\lib\\site-packages', 'C:\\OSGeo4W\\apps\\Python25\\lib\\site-packages\\win32', 'C:\\OSGeo4W\\apps\\Python25\\lib\\site-packages\\win32\\lib', 'C:\\OSGeo4W\\apps\\Python25\\lib\\site-packages\\Pythonwin', 'C:\\OSGeo4W\\apps\\Python25\\lib\\site-packages\\wx-2.8-msw-unicode', 'C:\\OSGeo4W\\apps\\qgis-dev\\python\\plugins\\fTools\\tools', 'E:/IMAGES_POOL/DTClass/Project']

[Voltron]

Общедоступная версия RasterCalc работает только с одноканальными изображениями. Т.е. в вашем случаем необходимо обработать каждый канал отдельно, а потом собрать их вместе при помощи gdal_merge

[SS_Rebelious]

А можно ли комбинировать несколько операций в одном выражении? Например, к пикселям со значением x прибавить значение y, а пиксели со значением > z присвоить значения другого растра? Или без промежуточных растров не обойтись?

UPD: использую например, такое выражение (не совсем такое, как в описании выше):
gt((([hydro_poly]@1=0) +1)@1, 10, 0)

Спойлер

Код: Выделить всё

An error has occured while executing Python code:

Traceback (most recent call last):
  File "/home/yury/.qgis//python/plugins/rastercalc/rastercalcdialog.py", line 269, in accept
    result = rastercalcengine.evaluateStack( tmpStack, iter, sizeX, blk )
  File "/home/yury/.qgis//python/plugins/rastercalc/rastercalcengine.py", line 209, in evaluateStack
    inRaster = evaluateStack( s, row, size, count )
  File "/home/yury/.qgis//python/plugins/rastercalc/rastercalcengine.py", line 223, in evaluateStack
    return returnBand( lay, num, row, size, count )
  File "/home/yury/.qgis//python/plugins/rastercalc/rastercalcengine.py", line 74, in returnBand
    return rasterUtils.getRasterBand( layerName, bandNum, row, size, count )
  File "/home/yury/.qgis//python/plugins/rastercalc/rastercalcutils.py", line 254, in getRasterBand
    gdalData = gdal.Open( str( rasterList[ name ].source() ) )
TypeError: unhashable type: 'numpy.ndarray'

Python version:
2.7.2 (default, Aug 19 2011, 20:41:43) [GCC]


QGIS version:
1.8.0-Lisboa Lisboa, exported

Python path: ['/home/yury/.qgis//python/plugins/sextante', '/home/yury/.qgis//python/plugins/sextante_taudem', '/usr/share/qgis/python', '/home/yury/.qgis//python', '/home/yury/.qgis//python/plugins', '/usr/share/qgis/python/plugins', '/windows/D/\xd0\x9c\xd0\xbe\xd0\xb8 \xd0\xb4\xd0\xbe\xd0\xba\xd1\x83\xd0\xbc\xd0\xb5\xd0\xbd\xd1\x82\xd1\x8b', '/usr/share/qgis/python/qgis', '/usr/share/qgis/python', '/usr/lib/python27.zip', '/usr/lib64/python2.7', '/usr/lib64/python2.7/plat-linux2', '/usr/lib64/python2.7/lib-tk', '/usr/lib64/python2.7/lib-old', '/usr/lib64/python2.7/lib-dynload', '/usr/lib64/python2.7/site-packages', '/usr/lib64/python2.7/site-packages/PIL', '/usr/local/lib64/python2.7/site-packages', '/usr/local/lib/python2.7/site-packages', '/usr/lib64/python2.7/site-packages/gtk-2.0', '/usr/lib/python2.7/site-packages', '/usr/lib/python2.7/site-packages/setuptools-0.6c11-py2.7.egg-info', '/home/yury/.qgis/python/plugins/imgboundary/logic', '/home/yury/.qgis/python/plugins/imgboundary/gui', '.', '/usr/share/qgis/python', '/usr/share/qgis/python/plugins/fTools/tools', '/home/yury/.qgis/python/plugins/openlayersov/logic', '/home/yury/.qgis/python/plugins/openlayersov/gui']

[Voltron]

А можно ли комбинировать несколько операций в одном выражении?

Можно. Например, расчет NDVI состоит из нескольких операций (вычитание, сложение, деление )

UPD: использую например, такое выражение (не совсем такое, как в описании выше):
gt((([hydro_poly]@1=0) +1)@1, 10, 0)

Тут что-то совсем страшное написано. В грамматике RasterCalc таких вещей как присваивание нет.