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

Tue Feb 10 20:48:56 EST 2009

Sorry, this was not the version I intended to send initially.
Here it is...
_______________________________________

# 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 = 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.dtype = numpy.int
	out.SetGeoTransform( geoT )
	out.SetProjection( proJ )
	driver.CreateCopy(output_filename, out)
	#SaveArray(result, output_filename, 'GTiff')
	#Empty memory objects
	del result, b, out

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

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 Rice Research Institute
Office: http://www.irri.org/gis
Perso: http://www.freewebs.com/ychemin
YiKingDo: http://yikingdo.unblog.fr/