On 01/31/2010 07:40 PM, Alexander Bruy wrote:
Hi all
I want to write a raster calculator and have two problems with Python
and GDAL.
My first problem: when calculating a difference between two raster bands I
get a wrong results. Raster in ERDAS IMAGINE format (HFA) and bands extracted
correctly. When I print each band and subtraction result I get something
similar to this
BAND 1
[[ 118 121 121 ..., 111 108 107]....]]
BAND 5
[[ 156 162 159 ..., 148 145 142]....]]
DIFFERENCE
[[218 215 218 ..., 219 219 221]....]]
But when calculate difference between two separate TIFFs (each band in
separate file) - result is ok
RASTER 1
[[ 118. 121. 121. ..., 111. 108. 107.]....]]
RASTER 5
[[ 156. 162. 159. ..., 148. 145. 142.]....]]
DIFFERENCE
[[-38. -41. -38. ..., -37. -37. -35.]....]]
Here is my code and sample files
http://www.4shared.com/file/212381308/a3e3f140/gdal_numpytar.html.
Run examples as:
- bands subtraction: test_bands.py landsat-6bands/s1986-06-12-163.img
- files subtraction: test_files.py landsat-separate-files/band1.tif
landsat-separate-files/band5.tif
Can anyone help? What is wrong in my code?
Hi Alexander,
the problem is that you're working with unsigned 8bit integers
(numpy.uint8 or numpy.byte). Therefore, the negative result (-38) is
interpreted as 256 - 38: 218. This is called 'integer overflow'
[http://en.wikipedia.org/wiki/Integer_overflow]. To allow for negative
values, you should do your calculations in a variable type that allows
for the storage of such values, e.g. in this case a signed 16bit integer
(numpy.int16) would suffice. For convenience sake, you could also do all
your calculations in floating point variables, so you won't be bitten by
integer overflows again. This will cost you more memory however, and in
some cases (not this one) you might loose some precision.
Vincent.
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
