Синтаксис растрового калькулятора gdal_calc для логических операторов и других функций

13

В документации по gdal_calc указан растровый калькулятор командной строки с пустым синтаксисом . Позже есть несколько примеров, где в одном из них:

gdal_calc.py -A input.tif --outfile = result.tif --calc = "A * (A> 0)" --NoDataValue = 0 - означает установить значения от нуля и ниже до нуля

К сожалению, нет примера для логических операторов, таких как:

--calc = "A * (A> 0 и A> B)" - означает оставить A, если A больше нуля и больше B, и установить остаток на ноль

Основываясь на логических функциях Numpy / Scipy, я ожидал бы написать логические операторы как:

--calc = "А * logical_and (А> 0, A> B)"

Я попробовал это, и это, кажется, работает, но я хотел бы быть уверен, что это правильно.

Аналогичным образом, если вы хотите минимум A и B:

--calc = "А * (А <= В) + В * (А> В)"

Вы можете просто написать:

--calc = "минимальный (А, В)"

Моя проблема в том, что я не могу найти кулинарную книгу, чтобы убедиться, что я правильно понял. Есть ли хорошая поваренная книга с продвинутыми примерами того, что можно и нельзя с gdal_calc?

Miro
источник

Ответы:

10

В источнике для gdal_calc.py, расчет сделан непосредственно, используя eval:

myResult = eval(opts.calc, global_namespace, local_namespace)

Это предполагает, что любое правильно сформированное выражение, которое также оценивается в командной строке, будет работать. Согласно документации, вы можете использовать gdalnumeric синтаксис с +-/*и / или numpyфункциями. Вы можете проверить свои функции с помощью небольших фиктивных массивов в интерактивной оболочке, а затем использовать те же вызовы в gdal_calc.

Имейте в виду, что объединение нескольких numpyфункций может создать временные массивы в памяти, которые могут существенно увеличить использование памяти, особенно при работе с большими изображениями.

Вы можете посмотреть в документации по списку всех функций: подпрограммы . Те, которые вы ищете, вероятно, здесь: математика или здесь: routines.logic .

Вот откуда берутся такие функции, как минимум, просто пространство имен уже импортировано. На самом деле, это numpy.minimum и т. Д.

Вениамин
источник
1
Спасибо, Бен, это еще один способ, о котором я понятия не имел. Еще после некоторой кулинарной книги, которая объяснила бы, что можно использовать, потому что eval не включает в себя функции минимума () и т. Д., Которые фактически можно использовать в выражении.
Миро
8

Исходя из ответа Бенджамина, вы можете использовать logic_or () или logic_and (). См. Http://docs.scipy.org/doc/numpy/reference/routines.logic.html . Следующий пример работал хорошо для меня. Это устанавливает все значения между 177 и 185 (включительно) равными 0, что затем обрабатывается как nodata.

gdal_calc.py -A input.tif --outfile=output.tif --calc="A*logical_or(A<=177,A>=185)" --NoDataValue=0
Tybion
источник
1

У меня был растр, где значения находились в диапазоне от -1 до 3, где ноль - действительное число. У меня были некоторые проблемы с созданием выражения gdal_calc, поэтому я сделал это быстрое и яростное решение.

#!/usr/bin/env python3

fileNameIn = "/tmp/geotiff/Global_taxonomic_richness_of_soil_fungi.tif"
fileNameOut = "/tmp/geotiff/Global_taxonomic_richness_of_soil_fungi.tiff"
dst_options = ['COMPRESS=DEFLATE',"PREDICTOR=3","TILED=YES"]
noDataValue = -3.4028234663852886e+38

from osgeo import gdal
import numpy

src_ds = gdal.Open(fileNameIn)
format = "GTiff"
driver = gdal.GetDriverByName(format)
dst_ds = driver.CreateCopy(fileNameOut, src_ds, False ,dst_options)

# Set location
dst_ds.SetGeoTransform(src_ds.GetGeoTransform())
# Set projection
dst_ds.SetProjection(src_ds.GetProjection())
srcband = src_ds.GetRasterBand(1)

dataraster = srcband.ReadAsArray().astype(numpy.float)
#Rplace the nan value with the predefiend noDataValue
dataraster[numpy.isnan(dataraster)]=noDataValue

dst_ds.GetRasterBand(1).WriteArray(dataraster)
dst_ds.GetRasterBand(1).SetNoDataValue(noDataValue)

dst_ds = None
Хорхе Мендес
источник