gdal_calc-Rasterrechnersyntax für logische Operatoren und andere Funktionen

13

In der Dokumentation für gdal_calc wird der Befehlszeilen-Rasterrechner mit Numpy-Syntax angegeben . Später gibt es einige Beispiele, wo in einem von ihnen:

gdal_calc.py -A input.tif --outfile = result.tif --calc = "A * (A> 0)" --NoDataValue = 0 - bedeutet , dass Werte von Null und darunter auf Null gesetzt werden

Leider gibt es kein Beispiel für logische Operatoren wie:

--calc = "A * (A> 0 und A> B)" - bedeutet dass A , wenn A größer als Null und größer als B ist, und der Rest auf Null gesetzt wird

Beyogen auf Numpy / Scipy-Logikfunktionen würde ich erwarten, logische Operatoren wie folgt zu schreiben:

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

Ich habe es versucht und es scheint zu funktionieren, aber ich möchte sicher sein, dass dies korrekt ist.

In ähnlicher Weise, wenn Sie ein Minimum von A und B wünschen:

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

Sie können einfach schreiben:

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

Mein Problem ist, dass ich kein Kochbuch finde, um sicherzustellen, dass ich das richtig mache. Gibt es ein gutes Kochbuch mit fortgeschrittenen Beispielen dafür, was mit gdal_calc möglich ist und was nicht?

Miro
quelle

Antworten:

10

In der Quelle für gdal_calc.py erfolgt die Berechnung direkt mit eval:

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

Das würde bedeuten, dass jeder wohlgeformte Ausdruck, der auch in der Befehlszeile ausgewertet wird, funktioniert. Gemäß der Dokumentation können Sie die gdalnumerische Syntax mit +-/*und / oder numpyFunktionen verwenden. Sie können Ihre Funktionen mit kleinen Dummy-Arrays in der interaktiven Shell testen und dann dieselben Aufrufe in gdal_calc verwenden.

Beachten Sie, dass die Verkettung mehrerer numpyFunktionen möglicherweise zu temporären In-Memory-Arrays führt, die die Speichernutzung erheblich erhöhen können, insbesondere bei großen Bildern.

In der numpy-Dokumentation finden Sie eine Liste aller Funktionen: Routinen . Die, nach denen Sie suchen, sind wahrscheinlich hier: math oder hier: routines.logic .

Hier kommen Funktionen wie Minimum her, nur dass der Namespace bereits importiert ist. Wirklich, es ist numpy.minimum usw

Benjamin
quelle
1
Danke Ben, das ist eine andere Möglichkeit, von der ich keine Ahnung hatte. Nach ein paar Kochbüchern, die erklären, was verwendet werden kann, weil eval keine minimum () etc-Funktionen enthält, die tatsächlich in expression verwendet werden können.
Miro
8

In Anlehnung an Benjamins Antwort können Sie logical_or () oder logical_and () verwenden. Siehe http://docs.scipy.org/doc/numpy/reference/routines.logic.html . Das folgende Beispiel hat für mich gut funktioniert. Dies setzt alle Werte zwischen 177 und 185 (einschließlich) auf 0, was dann als Knoten behandelt wird.

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

Ich hatte ein Raster mit Werten zwischen -1 und 3, wobei Null eine gültige Zahl ist. Ich hatte einige Probleme damit, einen gdal_calc-Ausdruck zu erstellen, so dass sich diese schnelle und wütende Lösung ergab.

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