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
