Hi,
I wrote a function to compute openness (a digital elevation model
attribute). The function opens the gdal-compatible input raster and
reads the values in a square window around each cell. The results are
stored in a "1 line buffer". At the end of the line, the buffer is
written to the output raster.
It runs very slowly and I suspect the raster access to be the main
bottleneck. Do you have any idea to speed-up the process ?
Thanks for any help.
Julien
Here is the code:
#################################################
## openness.py
##
## Julien Fiore, UNIGE
##
## code inspired by Matthew Perry's "slope.cpp"
# import
import numpy
import gdal
from gdal.gdalconst import * # for GA_ReadOnly in gdal.gdal.Open()
def openness(infile, outfile, R=1):
"""
Calculates the openness of a gdal-supported raster DEM.
Returns nothing.
Parameters:
infile: input file (georeferenced raster)
outfile: output file (GeoTIF)
R: opennness radius (in cells). Window size = (2*R +
1)**2
"""
# open inGrid
try:
inGrid = gdal.gdal.Open(infile, GA_ReadOnly )
except IOError:
print 'could not open inGrid'
inGridBand = inGrid.GetRasterBand(1) # get raster band
# get inGrid parameters
geoTrans = inGrid.GetGeoTransform() # returns list
(xOrigin,xcellsize,...)
cellSizeX = geoTrans[1]
cellSizeY = geoTrans[5]
nXSize = inGrid.RasterXSize
nYSize = inGrid.RasterYSize
nullValue = inGridBand.GetNoDataValue()
# Create openness output file
format = "GTiff"
driver = gdal.gdal.GetDriverByName( format )
outGrid=driver.CreateCopy(outfile,inGrid) # Creates output raster
with
# same properties as
input
outGridBand = outGrid.GetRasterBand(1) # Access Band 1
outGridBand.SetNoDataValue(nullValue)
# moving window
winSize= 2*R + 1
# openness buffer array (1 line)
buff = inGridBand.ReadAsArray(xoff=0, yoff=0, win_xsize=nXSize,
win_ysize=1)
# create distance array for Openness
dist=numpy.zeros((winSize,winSize),float)
for i in range(winSize):
for j in range(winSize):
dist[i][j]=numpy.sqrt((cellSizeX*(i-R))**2 +
(cellSizeY*(j-R))**2)
# openness loop
for i in range(nYSize):
for j in range(nXSize):
containsNull = 0
# excludes the edges
if (i <= (R-1) or j <= (R-1) or i >= (nYSize-R) or j >=
(nXSize-R) ):
buff[0][j] = nullValue
continue # continues with loop next iteration
# read window
win = inGridBand.ReadAsArray(j-R, i-R, winSize, winSize)
# Check if window has null value
for value in numpy.ravel(win):
if value == nullValue :
containsNull = 1
break # breaks out of the smallest enclosing loop
if (containsNull == 1):
# write Null value to buffer
buff[0][j] = nullValue
continue # continues with loop next iteration
else:
# openness calculation with "tan=opp/adj"
# East
E=180.0 # 180 = max angle possible between nadir and
horizon
for k in range(R+1,2*R):
angle= 90.0 -
numpy.arctan((win[k][R]-win[R][R])/dist[k][R])
if angle<E: E=angle
# West
W=180.0
for k in range(0,R-1):
angle= 90.0 -
numpy.arctan((win[k][R]-win[R][R])/dist[k][R])
if angle<W: W=angle
# North
N=180.0
for k in range(0,R-1):
angle= 90.0 -
numpy.arctan((win[R][k]-win[R][R])/dist[R][k])
if angle<N: N=angle
# South
S=180.0
for k in range(R+1,2*R):
angle= 90.0 -
numpy.arctan((win[R][k]-win[R][R])/dist[R][k])
if angle<S: S=angle
# mean openness angle
buff[0][j]= (E+W+N+S)/4
# Writes buffer to file
outGridBand.WriteArray(array=buff,xoff=0, yoff=i)
outGridBand.FlushCache()
# close datasets
del inGrid
del outGrid
return
if __name__ == "__main__":
# meazure time taken
from time import clock
startTime=clock() # START...
# Import Psyco if available
psycoBool=False
try:
import psyco
psyco.full()
psycoBool=True
except ImportError:
pass
startTime=clock()
# openness parameters
infile = 'E:/GIS/CH/VD/lidar/lowcut_filter/gaussian/xsmall5'
outfile = 'E:/GIS/CH/VD/lidar/openness/openness.tif' # geotiff
(.tif)
R=4
# call openness function
openness(infile,outfile,R)
endTime = clock() # ...END
# print information
print '\nOpenness OK, R=' +str(R) + '; Time elapsed: ' +
str(endTime-startTime) + ' seconds'
print 'psyco = ' + repr(psycoBool)
##############################################################
--
http://mail.python.org/mailman/listinfo/python-list
