How does gdal_calc numpy operators work?

Question

Gdal_calc manual

Gdal_calc is introduced as :

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

Unexpected results

I made some good uses of gdal_calc in the past, yet, recent results makes no sences to me. I have 3 input files hillshades_A.tmp.tif, hillshades_B.tmp.tif, hillshades_C.tmp.tif and one output composite.tif. I want to make calculations with these inputs. I will follow pixel x=10 y=10 with values are A=222, B=220, C=224 to check the correctness of the output. Let's start.

Example:

 # 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

Then the fun start !

# 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?

And I may go again and yet again. These are not simple basic arithmetic [...] such as +-*/ .

Question

How does gdal_calc operators works ? and... Is there a comprehensive manual for this gdal_calc / numppy syntaxes equations ? (I didn't find any solid one via google)

---

EDIT: may this behavior be the byproduct of [0-255] greyscale hillshade as input.

$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 seems to have more operators :

+ 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

I haven't found clear and proper examples for how the exotic operators should be used. If you have something, feel free to share.

Answer 1

As observed in the comments, what happens here is not due to the operators behaving strangely but to the fact that the output is written in Byte. So the computation is also done in Byte and every output is modulo 255. Even more, every temporary output is modulo 255, so you must take care of the precedence of the operator to make sure that no temporary result is greater than 255.

The simplest workaround is to convert your input data into integer with gdal_translate (output as VRT would not take a lot of space on your disk, but you can also write another tif using -ot Int16 ). And you should also define the output type to get an appropriate output type (that is, NOT Byte) for your operations.

--type=datatype output datatype, must be one of ['Int32', 'Int16', 'Float64', 'UInt16', 'Byte', 'UInt32', 'Float32']