sintaxe da calculadora raster gdal_calc para operadores lógicos e outras funções

Na documentação do gdal_calc , é apresentada a calculadora raster de linha de comando com sintaxe numpy . Mais tarde, existem poucos exemplos em que um deles:

gdal_calc.py -A input.tif --outfile = result.tif --calc = "A * (A> 0)" --NoDataValue = 0 - significa definir valores de zero e abaixo como nulos

Infelizmente, não há exemplo em operadores lógicos como:

--calc = "A * (A> 0 e A> B)" - significa manter A se A maior zero e maior B e definir o restante como nulo

Com base nas funções lógicas Numpy / Scipy, eu esperaria escrever operadores lógicos como:

--calc = "A * lógica_e (A> 0, A> B)"

Eu tentei isso e parece funcionar, mas gostaria de ter certeza de que está correto.

Da mesma forma, se você deseja um mínimo de A e B:

--calc = "A * (A <= B) + B * (A> B)"

Você pode simplesmente escrever:

--calc = "mínimo (A, B)"

Meu problema é que não consigo encontrar nenhum livro de receitas para garantir que esteja certo. Existe algum bom livro de receitas com exemplos avançados do que é e não é possível com o gdal_calc?

Na fonte de gdal_calc.py, o cálculo é feito diretamente usando eval:

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

Isso sugere que qualquer expressão bem formada que também avalie na linha de comando funcionará. De acordo com a documentação, você pode usar a sintaxe gdalnumeric com +-/*, e / ou numpyfunções. Você pode testar suas funções usando pequenas matrizes fictícias no shell interativo e, em seguida, use as mesmas chamadas em gdal_calc.

Lembre-se de que o encadeamento de várias numpyfunções provavelmente produzirá matrizes temporárias na memória que podem aumentar substancialmente o uso da memória, especialmente ao lidar com imagens grandes.

Você pode consultar a documentação numpy para obter uma lista de todas as funções: rotinas . Os que você procura provavelmente estão aqui: math ou aqui: routines.logic .

É daí que vêm as funções como mínimo, mas o espaço para nome já foi importado. Realmente, é numpy.minimum, etc

Seguindo a resposta de Benjamin, você pode usar logic_or () ou logic_and (). Consulte http://docs.scipy.org/doc/numpy/reference/routines.logic.html . O exemplo a seguir funcionou bem para mim. Isso define todos os valores entre 177 e 185 (inclusive) para 0, que é tratado como nodata.

gdal_calc.py -A input.tif --outfile=output.tif --calc="A*logical_or(A<=177,A>=185)" --NoDataValue=0

Eu tinha uma varredura em que os valores variavam entre -1 e 3, em que zero é um número válido. Eu tive alguns problemas ao criar uma expressão gdal_calc, então criei esta solução rápida e furiosa.

#!/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