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)

################################################## ############