[gdal-dev] Calculations with uint16 rasters

[Even Rouault]

On mercredi 18 septembre 2019 10:00:51 CEST Alexandre Neto wrote:
> Hi,
> 
> I am new to this list so sorry if this has ever been asked or discussed.
> 
> I am using two sentinel bands in gdal_calc. they are in jp2 format and
> uint16. I was surprised that subtracting one for the other (A - B) was
> giving me strange results.
> 
> gdal_calc.py --calc "(A-B)" --format GTiff --type Float32 -A
> /home/aneto/ARSI/data/S2A_MSIL2A_20170513T172301_N9999_R012_T14TPP_20190912T
> 155528.SAFE/GRANULE/L2A_T14TPP_A009876_20170513T172257/IMG_DATA/R10m/T14TPP_
> 20170513T172301_B08_10m.jp2 --A_band 1 -B
> /home/aneto/ARSI/data/S2A_MSIL2A_20170513T172301_N9999_R012_T14TPP_20190912T
> 155528.SAFE/GRANULE/L2A_T14TPP_A009876_20170513T172257/IMG_DATA/R10m/T14TPP_
> 20170513T172301_B04_10m.jp2 --B_band 1 --outfile
> /tmp/processing_90d61cb0e2a74784a9fcf4b11f9fd7dc/9945d3aec4a749d98f2a3d14694
> abe16/OUTPUT.tif
> 
> The result only returns positive integers, and the integers scale seems
> shifted to accommodate the negative values to something above 0.
> 
> Strangely the fix comes by multiplying the values by 0.1. So (A*1.0 -
> B*0.1) works as expected.
> 
> Is this a known issue? Is it the expected behavior?

It's kind of expected given how NumPy works, but a lot of people get into this 
trap.

Instead of multiplying by 1, you can also do A.astype(numpy.float32). This 
should be marginally faster.

This is very much the same as doing in C/C++

unsigned int A = 1;
unsigned int B = 2;
float C = A - B;

C will be 4294967296 (initialy I used 'unsigned short' for A and B, and I got 
-1 as the result, because in C/C++, when there's an arithmetic operation on 
types smaller than int, the operand get promoted to int first, which is 
signed)

We could possibly modify gdal_calc to grab the input raster on the specified 
output type, but I'm not completely sure this is a good idea, as there could 
be use cases where people would want the initial data type to be preserved.

But feel free to suggest a hint in the documentation of
https://gdal.org/programs/gdal_calc.html
I'd suggest both a note in the doc of the --type switch, and add an example 
similar to yours at the end.

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com