Sunday, 6 November 2016

gdal - How does gdal_calc numpy operators work?


Gdal_calc manual


Gdal_calc is introduced as :


gdal_calc.py - Command line raster calculator with numpy syntax
gdal_calc.py [-A ] [--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.




No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...