[gdal-dev] satellite image processing in Python

[Yann Chemin]

Hello,

Would need a small help to troubleshoot this small code:

1 - When the for loop is over, the last file is not completely closed
(~480Mb instead of 509Mb)
Is there a way to force it to close? (opening another file forces the
previous file to be closed with a proper"?" 509Mb)

2 -  A segmentation fault is crashing Python when using LoadFile() on
any one of the files generated in the loop.

Thank you,
Yann

-------------------------------------------------------------------------------------

# Landsat 7 ETM + image downlaoded from GLCF
#ftp.glcf.umiacs.umd.edu/glcf/Landsat/WRS2/p127/r049/p127r049_7x20001104.ETM-EarthSat-Orthorectified/
#p127r049 2000/11/04 around Ubon Ratchathani, North East Thailand.
# Metadata in the text file: p127r049_7x20001104.met

#Purpose: Calculate NDVI and Albedo
#Usage: /usr/bin/python -u "ndvi_albedo.py"

#if problem of libgrass: MAPSET is not set
#then: set GDAL_SKIP = GRASS

F=[]
F.append('p116r049_7t20010518_z51_nn10.tif')
F.append('p116r049_7t20010518_z51_nn20.tif')
F.append('p116r049_7t20010518_z51_nn30.tif')
F.append('p116r049_7t20010518_z51_nn40.tif')
F.append('p116r049_7t20010518_z51_nn50.tif')
F.append('p116r049_7t20010518_z51_nn70.tif')

# For Image Processing
from math import *
import numpy
from osgeo import gdalnumeric
from osgeo import gdal
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,241.1,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 = 51.3428710
doy = 311
#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):
	b = LoadFile(F[i])
	result = int(100*((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
	out = OpenArray(result)
	out.SetGeoTransform( geoT )
	out.SetProjection( proJ )
	driver.CreateCopy(output_filename, out)
	#Empty memory objects
	del result, b, out

#NDVI
NIR = LoadFile('b4_ref.tif') <------Seg Fault here
Red = LoadFile('b3_ref.tif')
ndvi = (NIR-Red)/(NIR+Red)


-- 
Yann Chemin
International Rice Research Institute
Office: http://www.irri.org/gis
Perso: http://www.freewebs.com/ychemin
YiKingDo: http://yikingdo.unblog.fr/