1

I was carrying out a calculation in raster calculator in ArcMap 10.4 where I multiplied two rasters together. Both rasters had the same extent, were snapped to each other and have the same pixels overlapping. When using the identify tool, I noticed that the result was different to if I worked it out myself. For example, the calculation for a particular pixel was 37.572238*90428390 (so should equal 3,397,596,991.03682). However the actual result from the raster calculator equaled 3,397,596,987.419. I also carried this out in QGIS to test and got another different result of 3,397,596,928.000.

Does anyone know why this happens? I thought it may have something to do with the precision as the numbers are so large...I know it isn't a huge difference but I need the result to be as accurate as possible.

0

You're right, this is likely to be a floating point precision issue. Floating point numbers are inherently 'lossy' as you're trying to represent an infinite number space in a fixed number of bits.

This will depend on

  • your operating system and software (are both 32 bit or 64 bit)
  • the data type in the raster file (32 bit float or 64 bit float?). You can use the gdalinfo command line tool to find out.

for example, in Python you can see different results by forcing it to use 32 or 64 bit floats

>>> import numpy as np
>>> np.float32(37.572238)*np.int(90428390)
3397597074.4225693
>>> np.float64(37.572238)*np.int(90428390)
3397596991.0368199

When I try vanilla python or the Mac calculator, I get the second version

Your QGIS version that returns an 'integer float' (3,397,596,928.000) may be the result of casting the result to an integer if one raster is float and the other is integer, then casting back to a float. That doesn't sound correct but it might explain that odd result.

You can get around this one of two ways

  • use a 64 bit TIFF
  • use 64bit integer, multiply your floats by 1e6, then divide the result by 1e6. This is rounded to 4dp.
  • Thanks for your reply, I used both 64bit TIFFs but in both ArcMap and QGIS the output is still different, and the raster produced is set as a 32bit float data type- I don't think there is an option to set the desired data type, do you think it would be better to do this in python? – Issy Sep 13 '18 at 9:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.