gdal_calc sintaxis de la calculadora raster para operadores lógicos y otras funciones

Question

En documentación de gdal_calc se afirma Calculadora de rasterización en línea de comandos con sintaxis numpy . Más adelante hay algunos ejemplos en los que en uno de ellos:

gdal_calc.py -A input.tif --outfile=result.tif --calc="A*(A>0)" --NoDataValue=0 - significa establecer valores de cero e inferiores a null

Lamentablemente no hay ningún ejemplo sobre operadores lógicos como:

--calc="A*(A>0 y A>B)"- significa mantener A si A es mayor que cero y mayor que B y establecer el resto como nulo

Basado en Funciones lógicas Numpy/Scipy Yo esperaría escribir operadores lógicos como:

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

He probado esto y parece que funciona pero me gustaría tener la seguridad de que es correcto.

De la misma manera si quieres el mínimo de A y B:

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

Puedes simplemente escribir:

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

Mi problema es que no encuentro ningún libro de cocina para asegurarme de que lo hago bien. ¿Existe algún buen libro de recetas con ejemplos avanzados de lo que es y no es posible con gdal_calc?

Answer 1

En el código fuente de gdal_calc.py, el cálculo se realiza directamente utilizando eval :

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

Eso sugeriría que cualquier expresión bien formada que también se evalúe en la línea de comandos funcionará. Según la documentación, puede utilizar la sintaxis gdalnumeric con +-/* y/o numpy funciones. Puede probar sus funciones utilizando pequeñas matrices ficticias en el intérprete de comandos interactivo y, a continuación, utilizar las mismas llamadas en gdal_calc.

Tenga en cuenta que encadenar varios numpy es probable que produzca matrices temporales en memoria que pueden aumentar sustancialmente el uso de memoria, especialmente cuando se trata de imágenes de gran tamaño.

Puedes consultar la documentación de numpy para obtener una lista de todas las funciones: rutinas . Los que busca probablemente estén aquí: matemáticas o aquí: rutinas.lógicas .

De aquí salen funciones como mínimo, sólo que el espacio de nombres ya está importado. Realmente, es numpy.minimum, etc

Answer 2

Siguiendo con la respuesta de Benjamin, puedes utilizar logical_or() o logical_and(). Véase http://docs.scipy.org/doc/numpy/reference/routines.logic.html . El siguiente ejemplo me ha funcionado muy bien. Esto establece todos los valores entre 177 y 185 (inclusive) a 0, que luego se trata como nodata.

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

Answer 3

Tenía una trama en la que los valores oscilaban entre -1 y 3, donde cero es un número válido. Tuve algunos problemas haciendo una expresión gdal_calc así que hice esta solución rápida y 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

Answer 4

Esta es otra forma que funciona para hacerlo. Estoy usando GDAL vía QGIS en PyQGIS. Se puede especificar cualquier cosa en la 'FORMULA'.

surfaceA = 'P:/70210 Severn Trent Groundwater Map/GIS/OS_100k_WaterFeatures_Tamworth_DTM_interpolated_50mRes_WOCanal.tif'
surfaceB = 'P:/70210 Severn Trent Groundwater Map/GIS/Main_river_points_DTM_inTamworth_interpolated_50mRes.tif',
processing.run("gdal:rastercalculator", 
    {'INPUT_A': surfaceA,
     'BAND_A':1,
     'INPUT_B': surfaceB, 
     'BAND_B':1,
     'FORMULA':'numpy.min((A,B),axis=0)',
     'NO_DATA':None,
     'RTYPE':5,
     'OPTIONS':'',
     'EXTRA':'',
     'OUTPUT':'TEMPORARY_OUTPUT'
     })