Use of Python with GDAL. How to speed up ?

[Julien Fiore]

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)

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