[gdal-dev] Re: satellite image processing in Python

[Yann Chemin]

Hello,

Found out answer to question 1 it seems.

By naming the object for each action, then deleting it at the end of
the loop, it seems that the file is being written properly.

	geot=result.SetGeoTransform( geoT )
	proj=result.SetProjection( proJ )
	output=driver.CreateCopy(output_filename, result)
        del output, geot, proj

still loading the created file makes a segfault in ReadArray it seems.
I am using 1.5.2 for Lenny. Maybe should try 1.6

Yann
----------------------------------------------------------------
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'
	else:
		output_filename='b'+str(i+2)+'_ref'
	#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
-------------------------

2009/2/11 Yann Chemin <yann.chemin at gmail.com>:
> 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/
>



-- 
Yann Chemin
International Center of Water for Food Security
Office: http://www.icwater.org
Perso: http://www.freewebs.com/ychemin
YiKingDo: http://yikingdo.unblog.fr/