[gdal-dev] Raster algebra with Gdal and Python

[Vincent Schut]

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.