<html dir="ltr"><head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<style title="owaParaStyle"><!--P {
MARGIN-TOP: 0px; MARGIN-BOTTOM: 0px
}
--></style>
</head>
<body ocsi="x">
<div dir="ltr"><font face="Tahoma" color="#000000" size="2">Hello,</font></div>
<div dir="ltr"><font face="tahoma" size="2">I would like to use this code written by Yann Chemin (see Feb 2009 gdal-dev archives)</font></div>
<div dir="ltr"><font face="tahoma" size="2">but I keep getting a memory error. Can anyone help me fix this. I have put the entire error below the code.</font></div>
<div dir="ltr"><font face="tahoma" size="2">Thanks in advance</font></div>
<div dir="ltr"><font face="tahoma" size="2"></font> </div>
<div dir="ltr"><font face="Tahoma" color="#000000" size="2">F=[]<br>
F.append('e:\Imagery\Landsat\Landsat7\L71040036_03620071224_B10.TIF')<br>
F.append('e:\Imagery\Landsat\Landsat7\L71040036_03620071224_B20.TIF')<br>
F.append('e:\Imagery\Landsat\Landsat7\L71040036_03620071224_B30.TIF')<br>
F.append('e:\Imagery\Landsat\Landsat7\L71040036_03620071224_B40.TIF')<br>
F.append('e:\Imagery\Landsat\Landsat7\L71040036_03620071224_B50.TIF')<br>
F.append('e:\Imagery\Landsat\Landsat7\L72040036_03620071224_B70.TIF')</font></div>
<div dir="ltr"><font face="Tahoma" color="#000000" size="2"># For Image Processing<br>
from math import *<br>
import numpy<br>
from osgeo import gdal<br>
from osgeo import gdalnumeric<br>
from osgeo.gdal_array import *<br>
from osgeo.gdalconst import *</font></div>
<div dir="ltr"><font face="Tahoma" color="#000000" size="2"># Set our output file format driver<br>
driver = gdal.GetDriverByName( 'GTiff' )</font></div>
<div dir="ltr"><font face="Tahoma" color="#000000" size="2"># Set our output files Projection parameters<br>
# from input file number 1<br>
tmp = gdal.Open(F[0])<br>
geoT = tmp.GetGeoTransform()<br>
proJ = tmp.GetProjection()<br>
del tmp</font></div>
<div dir="ltr"><font face="Tahoma" color="#000000" size="2">#DN to radiance conversion factors<br>
#(Check if they are correct for your image)<br>
LMax = [191.6,196.5,152.9,157.4,31.060,10.8]<br>
Lmin = [-6.2,-6.4,-5.0,-5.1,-1.0,-0.35]<br>
#Exo-Atmospheric band-wise irradiance<br>
KEXO = [1970.0,1843.0,1555.0,1047.0,227.1,80.53]<br>
#Unless you use p127r049 2000/11/04, these must be changed<br>
sun_elevation = 27.7636885<br>
doy = 358<br>
#Cos(sun zenith angle) / (pi * (Sun-Earth distance)^2)<br>
constant=(cos((90-sun_elevation) * pi/180 ) / (pi * (1+0.01672 * sin<br>
(2 * pi * (doy - 93.5) / 365)) ** 2))</font></div>
<div dir="ltr"><font face="Tahoma" color="#000000" size="2">#Top Of Atmosphere Reflectance<br>
for i in range (0,6):<br>
tmp = gdal.Open(F[i])<br>
btmp = tmp.GetRasterBand(1)<br>
b = BandReadAsArray(btmp,0,0,tmp.RasterXSize,tmp.RasterYSize)<br>
result = (b * ((LMax[i]-Lmin[i])/255.0) + Lmin[i]) / (constant * <br>
KEXO[i])<br>
if i != 5:<br>
output_filename='b'+str(i+1)+'_ref.tif'<br>
else:<br>
output_filename='b'+str(i+2)+'_ref.tif'<br>
#Create output file with projection<br>
result = OpenArray(result)<br>
geot=result.SetGeoTransform( geoT )<br>
proj=result.SetProjection( proJ )<br>
output=driver.CreateCopy(output_filename, result)<br>
#Empty memory objects<br>
del result, b, tmp, btmp, output, geot, proj</font></div>
<div dir="ltr"><font face="Tahoma" color="#000000" size="2"></font> </div>
<div dir="ltr"><font face="Tahoma" color="#000000" size="2">#NDVI<br>
NIR = LoadFile('b4_ref.tif')<br>
Red = LoadFile('b3_ref.tif')<br>
ndvi = (NIR-Red)/1.0*(NIR+Red)</font></div>
<div dir="ltr"><font face="Tahoma" color="#000000" size="2"></font> </div>
<div dir="ltr"><font face="times new roman"></font> </div>
<div dir="ltr">Traceback (most recent call last):<br>
File "E:\python processing\Landsat_processing.py", line 75, in <module><br>
Red = LoadFile('b3_ref.tif')<br>
File "C:\Python25\lib\site-packages\osgeo\gdal_array.py", line 73, in LoadFile<br>
return DatasetReadAsArray( ds, xoff, yoff, xsize, ysize )<br>
File "C:\Python25\lib\site-packages\osgeo\gdal_array.py", line 90, in DatasetReadAsArray<br>
return BandReadAsArray( ds.GetRasterBand(1), xoff, yoff, xsize, ysize)<br>
File "C:\Python25\lib\site-packages\osgeo\gdal_array.py", line 138, in BandReadAsArray<br>
buf_xsize, buf_ysize, datatype )<br>
File "C:\Python25\Lib\site-packages\osgeo\gdal.py", line 760, in ReadRaster<br>
return _gdal.Band_ReadRaster(*args, **kwargs)<br>
MemoryError</div>
</body>
</html>