hi Even, If you think it's better to keep it that way, let's leave it as it is. :-)
I will definitely try to add that as a hint to the docs. Not sure if that's how you guys work, but I have created an issue for that taks: https://github.com/OSGeo/gdal/issues/1848 Thanks for the fast answer. Alex On Wed, Sep 18, 2019 at 10:19 AM Even Rouault <[email protected]> wrote: > 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 [email protected] https://lists.osgeo.org/mailman/listinfo/gdal-dev
