[gdal-dev] satellite image processing in Python

[Harris, Sarah L]

Hello,
I would like to use this code written by Yann Chemin (see Feb 2009 gdal-dev archives)
but I keep getting a memory error. Can anyone help me fix this. I have put the entire error below the code.
Thanks in advance

F=[]
F.append('e:\Imagery\Landsat\Landsat7\L71040036_03620071224_B10.TIF')
F.append('e:\Imagery\Landsat\Landsat7\L71040036_03620071224_B20.TIF')
F.append('e:\Imagery\Landsat\Landsat7\L71040036_03620071224_B30.TIF')
F.append('e:\Imagery\Landsat\Landsat7\L71040036_03620071224_B40.TIF')
F.append('e:\Imagery\Landsat\Landsat7\L71040036_03620071224_B50.TIF')
F.append('e:\Imagery\Landsat\Landsat7\L72040036_03620071224_B70.TIF')
# For Image Processing
from math import *
import numpy
from osgeo import gdal
from osgeo import gdalnumeric
from osgeo.gdal_array import *
from osgeo.gdalconst import *
# Set our output file format driver
driver = gdal.GetDriverByName( 'GTiff' )
# Set our output files Projection parameters
# from input file number 1
tmp = gdal.Open(F[0])
geoT = tmp.GetGeoTransform()
proJ = tmp.GetProjection()
del tmp
#DN to radiance conversion factors
#(Check if they are correct for your image)
LMax = [191.6,196.5,152.9,157.4,31.060,10.8]
Lmin = [-6.2,-6.4,-5.0,-5.1,-1.0,-0.35]
#Exo-Atmospheric band-wise irradiance
KEXO = [1970.0,1843.0,1555.0,1047.0,227.1,80.53]
#Unless you use p127r049 2000/11/04, these must be changed
sun_elevation = 27.7636885
doy = 358
#Cos(sun zenith angle) / (pi * (Sun-Earth distance)^2)
constant=(cos((90-sun_elevation) * pi/180 ) / (pi * (1+0.01672 * sin
(2 * pi * (doy - 93.5) / 365)) ** 2))
#Top Of Atmosphere Reflectance
for i in range (0,6):
        tmp = gdal.Open(F[i])
        btmp = tmp.GetRasterBand(1)
        b = BandReadAsArray(btmp,0,0,tmp.RasterXSize,tmp.RasterYSize)
        result = (b * ((LMax[i]-Lmin[i])/255.0) + Lmin[i]) / (constant *
KEXO[i])
        if i != 5:
                output_filename='b'+str(i+1)+'_ref.tif'
        else:
                output_filename='b'+str(i+2)+'_ref.tif'
        #Create output file with projection
        result = OpenArray(result)
        geot=result.SetGeoTransform( geoT )
        proj=result.SetProjection( proJ )
        output=driver.CreateCopy(output_filename, result)
        #Empty memory objects
        del result, b, tmp, btmp, output, geot, proj

#NDVI
NIR = LoadFile('b4_ref.tif')
Red = LoadFile('b3_ref.tif')
ndvi = (NIR-Red)/1.0*(NIR+Red)


Traceback (most recent call last):
  File "E:\python processing\Landsat_processing.py", line 75, in <module>
    Red = LoadFile('b3_ref.tif')
  File "C:\Python25\lib\site-packages\osgeo\gdal_array.py", line 73, in LoadFile
    return DatasetReadAsArray( ds, xoff, yoff, xsize, ysize )
  File "C:\Python25\lib\site-packages\osgeo\gdal_array.py", line 90, in DatasetReadAsArray
    return BandReadAsArray( ds.GetRasterBand(1), xoff, yoff, xsize, ysize)
  File "C:\Python25\lib\site-packages\osgeo\gdal_array.py", line 138, in BandReadAsArray
    buf_xsize, buf_ysize, datatype )
  File "C:\Python25\Lib\site-packages\osgeo\gdal.py", line 760, in ReadRaster
    return _gdal.Band_ReadRaster(*args, **kwargs)
MemoryError
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20090414/b5a47b69/attachment.html