[gdal-dev] Calculations with uint16 rasters

[Pedro Venâncio]

Hi Even,

> 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.

I think this should be considered at some time.

Of course this can break some current workflows, but that can be overcome,
adding a wildcard to --type=<datatype> to 'Keep the input datatype', and
keeping that option as the default.

Python itself already changed that behaviour in v3:

Python 2.7.14 (v2.7.14:84471935ed, Sep 16 2017, 20:25:58) [MSC v.1500 64
bit (AMD64)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> A=1
>>> B=2
>>> A-B
-1
>>> A/B
0

Python 3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit
(AMD64)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> A=1
>>> B=2
>>> A-B
-1
>>> A/B
0.5

and so, I believe it can be done with gdal_calc.

If a user choose a --type=Float32, it's almost sure that he wants the
result as a decimal with a considerable level of precision, and the cast of
the inputs should be done. If the user want the initial data type be
preserved, he could simply omit the --type parameter, and 'Keep the input
datatype' option was applied.

I think this would be more coherent and more understandable for users less
used to computer languages.

Thank you very much!

Best regards,
Pedro Venâncio

Even Rouault <even.rouault at spatialys.com> escreveu no dia quarta,
18/09/2019 à(s) 10:19:

> 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
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20190918/1ab8135e/attachment.html>