¿Cómo funcionan los operadores numpy de gdal_calc?

Question

Manual de Gdal_calc

Gdal_calc se presenta como :

gdal_calc.py - Command line raster calculator with numpy syntax
gdal_calc.py [-A <filename>] [--A_band] [-B...-Z filename] [other_options]
OPTIONS:
--calc=CALC           calculation in gdalnumeric syntax using +-/* or any
                      numpy array functions (i.e. logical_and())
DESCRIPTION
Command line raster calculator with numpy syntax. Use any basic arithmetic
supported by numpy arrays such as +-*/ along with logical operators such as >.
EXAMPLE
add two files together
   gdal_calc.py -A input1.tif -B input2.tif --outfile=result.tif --calc="A+B"
average of two layers
   gdal_calc.py -A input.tif -B input2.tif --outfile=result.tif --calc="(A+B)/2"`
set values of zero and below to null
   gdal_calc.py -A input.tif --outfile=result.tif --calc="A*(A>0)" --NoDataValue=0

Resultados inesperados

He hecho algunos buenos usos de gdal_calc en el pasado, sin embargo, los resultados recientes no tienen sentido para mí. Tengo 3 archivos de entrada hillshades_A.tmp.tif , hillshades_B.tmp.tif , hillshades_C.tmp.tif y una salida composite.tif . Quiero hacer cálculos con estas entradas. Voy a seguir pixel x=10 y=10 con valores son A= 222 , B= 220 , C= 224 para comprobar que la salida es correcta. Empecemos.

Ejemplo:

 # ASSIGN VALUE. Expect: 300.
 $gdal_calc.py -A ./hillshades_A.tmp.tif --outfile=./composite.tif  \
   --calc="300" --NoDataValue=-1
 $gdallocationinfo ./composite.tif 10 10 -valonly
 >255          # Ergh? Ok, as it's a grey scale [0-255], such auto-correction is ok

Entonces empieza la diversión.

# SUM: A+B or (A+B), when A=222, B=220. Expect: 442 or 255.
gdal_calc.py -A hillshades_A.tmp.tif -B hillshades_B.tmp.tif --outfile=./composite.tif \
   --calc="(A+B)" --NoDataValue=-1
$gdallocationinfo ./composite.tif 10 10 -valonly
>186          # Ergh? This is A+B-256=442-256=186

# AVERAGE: (A+B)/2, when A=222, B=220. Expect: 221 or 255.
gdal_calc.py -A hillshades_A.tmp.tif -B hillshades_B.tmp.tif --outfile=./composite.tif  \
   --calc="(A+B)/2" --NoDataValue=-1
$gdallocationinfo ./composite.tif 10 10 -valonly
>93           # Ergh? ok, per previous 186/2=93. 
# must use (A/2)+(B/2) => 221        

# MULTIPLY: A*B, when A=222, B=220. Expect: 48840 or 255
gdal_calc.py -A hillshades_A.tmp.tif -B hillshades_B.tmp.tif --outfile=./composite.tif  \
   --calc="A*B" --NoDataValue=-1
$gdallocationinfo ./composite.tif 10 10 -valonly
>200           # Ergh? ok, it's modulo: 222*220 % 256=200  

# MULTIPLY+SCALE BACK: A*B/255, when A=222, B=220. Expect: 191.5 rounded to 191 or 192.
gdal_calc.py -A hillshades_A.tmp.tif -B hillshades_B.tmp.tif --outfile=./composite.tif  \
   --calc="A*B/255" --NoDataValue=-1
$gdallocationinfo ./composite.tif 10 10 -valonly
>0             # Ergh? Maybe A*B/255=200/255 (per previous point) being <0 it get rounded to 0 firt. It confirms: Ergh?

Y puede que vaya otra vez y otra más. No son simples *_aritmética básica [...] como `+-/` ._**

Pregunta

¿Cómo funcionan los operadores de gdal_calc? y... ¿Existe un manual completo para estas ecuaciones de gdal_calc / numppy syntaxes? (No he encontrado ninguna sólida a través de Google)

---

EDITAR: puede ser este comportamiento el subproducto de [0-255] sombra de colina en escala de grises como entrada.

$gdalinfo ./hillshades_A.tmp.tif -mm
Driver: GTiff/GeoTIFF
Files: ./hillshades_A.tmp.tif
Size is 1920, 1950
Coordinate System is `'
Origin = (66.991666666666674,37.508333333333354)
Pixel Size = (0.016666666666667,-0.016666666666667)
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  66.9916667,  37.5083333) 
Lower Left  (  66.9916667,   5.0083333) 
Upper Right (  98.9916667,  37.5083333) 
Lower Right (  98.9916667,   5.0083333) 
Center      (  82.9916667,  21.2583333) 
Band 1 Block=1920x4 Type=Byte, ColorInterp=Gray
    Computed Min/Max=1.000,255.000
  NoData Value=0

---

Gdal_calc 's numpy parece tener más operadores :

+ addition 
- subtraction 
/ division 
* multiplication 
= equals to 
< less than 
> larger than 
! not equal to 
? if clause 
M maximum of two values 
m minimum of two values 
B bit level operator

No he encontrado ejemplos claros y adecuados de cómo se deben utilizar los operadores exóticos. Si tienes algo, no dudes en compartirlo.

Answer 1

Como se observa en los comentarios, lo que ocurre aquí no se debe a que los operadores se comporten de forma extraña, sino a que la salida está escrita en Byte. Así que el cálculo también se hace en Byte y cada salida es módulo 255. Es más, cada salida temporal es módulo 255, por lo que hay que cuidar la precedencia del operador para asegurarse de que ningún resultado temporal es mayor que 255.

La solución más sencilla es convertir los datos de entrada en enteros con gdal_translate (la salida como VRT no ocuparía mucho espacio en tu disco, pero también puedes escribir otro tif usando -ot Int16 ). Y también deberías definir el tipo de salida para obtener un tipo de salida apropiado (es decir, NO Byte) para tus operaciones.

--type=tipo de datos de salida, debe ser uno de ['Int32', 'Int16', 'Float64', 'UInt16', 'Byte', 'UInt32', 'Float32']