Syntaxe de la calculatrice raster gdal_calc pour les opérateurs logiques et autres fonctions

13

Dans la documentation de gdal_calc, il est indiqué Calculatrice raster en ligne de commande avec syntaxe numpy . Plus tard, il y a quelques exemples où dans l'un d'eux:

gdal_calc.py -A input.tif --outfile = result.tif --calc = "A * (A> 0)" --NoDataValue = 0 - signifie que les valeurs définies de zéro et inférieures sont nulles

Malheureusement, il n'y a pas d'exemple sur les opérateurs logiques comme:

--calc = "A * (A> 0 et A> B)" - signifie garder A si A plus grand zéro et plus grand B et mettre le reste à null

Sur la base des fonctions logiques Numpy / Scipy, je m'attendrais à écrire des opérateurs logiques comme:

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

J'ai essayé et cela semble fonctionner mais je voudrais être sûr que c'est correct.

De la même manière si vous voulez un minimum de A et B:

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

Vous pouvez simplement écrire:

--calc = "minimum (A, B)"

Mon problème est que je ne trouve aucun livre de cuisine pour m'assurer de bien faire les choses. Existe-t-il un bon livre de cuisine avec des exemples avancés de ce qui est et n'est pas possible avec gdal_calc?

Miro
la source

Réponses:

10

Dans la source de gdal_calc.py, le calcul se fait directement en utilisant eval:

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

Cela suggère que toute expression bien formée qui évalue également sur la ligne de commande fonctionnera. Selon la documentation, vous pouvez utiliser la syntaxe gdalnumérique avec +-/*et / ou les numpyfonctions. Vous pouvez tester vos fonctions à l'aide de petits tableaux factices dans le shell interactif, puis utiliser les mêmes appels dans gdal_calc.

Gardez à l'esprit que l'enchaînement de plusieurs numpyfonctions est susceptible de produire des matrices temporaires en mémoire qui peuvent augmenter considérablement l'utilisation de la mémoire, en particulier lorsqu'il s'agit de grandes images.

Vous pouvez consulter la documentation numpy pour une liste de toutes les fonctions: routines . Ceux que vous recherchez sont probablement ici: mathématiques ou ici: routines.logic .

C'est de là que viennent les fonctions comme minimum, c'est juste que l'espace de noms est déjà importé. Vraiment, c'est numpy.minimum, etc.

Benjamin
la source
1
Merci Ben, c'est une autre façon dont je n'avais aucune idée. Toujours après un livre de cuisine qui expliquerait ce qu'il est possible d'utiliser, car eval n'inclut pas les fonctions minimum (), etc., qui sont réellement possibles à utiliser dans l'expression.
Miro
8

À la suite de la réponse de Benjamin, vous pouvez utiliser logic_or () ou logic_and (). Voir http://docs.scipy.org/doc/numpy/reference/routines.logic.html . L'exemple suivant a très bien fonctionné pour moi. Cela définit toutes les valeurs comprises entre 177 et 185 (inclus) sur 0, qui est ensuite traitée comme nodata.

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

J'avais un raster où les valeurs variaient entre -1 et 3 où zéro est un nombre valide. J'ai eu quelques problèmes pour créer une expression gdal_calc, alors j'ai fait cette solution rapide et furieuse.

#!/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
Jorge Mendes
la source