# [gdal-dev] Raster algebra with Gdal and Python

Alexander Bruy alexander.bruy at gmail.com
Sun Jan 31 13:40:22 EST 2010

```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?

Second problem: processing large rasters. I get error

File "C:/Documents and Settings/alex/.qgis//python/plugins\rastercalc\rastercalcutils.py", line 162, in layerAsArray
File "C:\OSGeo4W\apps\gdal-16\pymod\osgeo\gdal.py", line 603, in ReadAsArray
return gdalnumeric.DatasetReadAsArray( self, xoff, yoff, xsize, ysize )
File "C:\OSGeo4W\apps\gdal-16\pymod\osgeo\gdal_array.py", line 108, in DatasetReadAsArray
return numpy.concatenate( array_list )
MemoryError

Large rasters is landsats with 6 bands, 30m/pixel, scene size is aprox. 185x185 km, aprox. 9000x9000 pixels. On disk one file occupied ~500-700 Mb. Some files
can be larger, up to 15000x15000 pix. and more.
Is it possible to work with large files? Maybe it is possible to process
small pieces of this files and then merge this pieces into one large file?