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[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)

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

 

[Caleb Hattingh]

Hi

The documentation for the python profiler in the python library
reference is extremely readable (well done James Roskin!?).

Profile your code, and when you find where the speed problem occurs,
try pitching just that section of code in comp.lang.python. You will
likely get much feedback. Everyone loves small optimization exercises,
and there be bots here who get the optimal solution on the first try.

Though I would like to (seems an interesting application), I don't have
time available to go through all your code now.

Thanks
Caleb

 

[Serge Orlov]

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 ?

Psyco doesn't understand numpy AFAIK. So you have you have two options:
1) Rewrite the code to use array.array from stdlib instead of numpy,
then (hopefully) psyco will optimize it. Since array.array is
one-dimensional, it's not going to be fun.
2) Rewrite the code to be vectorized (don't use psyco) Right now your
code *doesn't* get any speed benefit from numpy

-- Serge

 

[Julien Fiore]

Thanks Caleb for the advice,

I profiled my code with the following lines:

import profile, pstats

profile.runctx('openness(infile,outfile,R)',globals(),locals(),'profile.log')
p = pstats.Stats('profile.log')
p.sort_stats('time')
p.print_stats(10)

The outputs tells me that the openness() function takes most of the
time. This is not very useful for me, as openness() contains most of my
code. How to know which instructions take most time in the openness
function ? Do I have to separate the code into smaller functions ?

 

[Julien Fiore]

Thanks for the psyco information, Serge.

2) Rewrite the code to be vectorized (don't use psyco) Right now your
code *doesn't* get any speed benefit from numpy

I do not understand this point. How to rewrite the code ? Do you mean
in C ?

Julien

 

[Kent Johnson]

Thanks Caleb for the advice,

I profiled my code with the following lines:

import profile, pstats

profile.runctx('openness(infile,outfile,R)',globals(),locals(),'profile.log')
p = pstats.Stats('profile.log')
p.sort_stats('time')
p.print_stats(10)

The outputs tells me that the openness() function takes most of the
time. This is not very useful for me, as openness() contains most of my
code. How to know which instructions take most time in the openness
function ? Do I have to separate the code into smaller functions ?

Yes. The profiler collects statistics for functions, not for individual
lines of code.

Kent

 

[Serge Orlov]

Thanks for the psyco information, Serge.


I do not understand this point. How to rewrite the code ?

vectorized code means one operation work on multiple items. Here is
three variants of the most easy part of your function (i'm sure this
part is not hot, it's just an example, the hottest part I'm pretty sure
is the *whole* openness loop)

setup = """
import numpy, array
R = 4
winSize= 2*R + 1
cellSizeX = 11
cellSizeY = 23
"""
code1 = """
dist=numpy.zeros((winSize,winSize),float)
for i in range(winSize):
for j in range(winSize):
dist[j]=numpy.sqrt((cellSizeX*(i-R))**2
+(cellSizeY*(j-R))**2)
"""
code2 = """
ind = numpy.indices((winSize,winSize))
dist = numpy.sqrt((cellSizeX*(ind[0]-R))**2 +(cellSizeY*(ind[1]-R))**2)
"""
code3 = """
dist=array.array('f',chr(0)*(4*winSize*winSize))
for i in range(winSize):
for j in range(winSize):
dist[i*winSize+j]=((cellSizeX*(i-R))**2
+(cellSizeY*(j-R))**2)**0.5
"""

The first one is original, the second is vectorized, the last one is
using array.array. Here is the timing in microseconds for R=1,2,5,6,50

68.2184467579
95.1533784322
17.6263407779

173.079790867
108.401514718
39.155870369
802.235479649
174.977043172
161.933094849
1126.00120965
206.326781756
221.849145632
67636.3336681
6736.00415695
12675.9099271

As you see if R >=6 vectorized version is the fastest, otherwise
array.array is faster. You should also try psyco + array.array

Do you mean in C ?

I didn't mean it, but it's also an option if you need the best speed.
You can try pyrex first, before resorting to raw C.

Serge.

 

[Julien Fiore]

Thank you Serge for this generous reply,

Vectorized code seems a great choice to compute the distance. If I
understand well, vectorized code can only work when you don't need to
access the values of the array, but only need to know the indices. This
works well for the distance, but I need to access the array values in
the openness function.

As regards array.array, it seems a bit complicated to reduce all my 2D
arrays to 1D arrays because I am working with the gdal library, which
works only with 'numeric' arrays.

Julien

 

[Serge Orlov]

Thank you Serge for this generous reply,

Vectorized code seems a great choice to compute the distance. If I
understand well, vectorized code can only work when you don't need to
access the values of the array, but only need to know the indices.

You can access array elements in any order if you wrap your head around
how indexing array with arrays works. Here is I'm transposing matrix on
the fly

This
works well for the distance, but I need to access the array values in
the openness function.

What really kills vectorization is conditional control flow
constructions (if, break, continue). Some of them can be still be
converted if you change algorithm. But some kill vectorization. For
example first continue where you exclude edges is avoidable: you need
to break you main loop in two, first loop will initialize edges, the
second loop will work on the main image. But the second continue is a
killer. Another example of algorightm change is small loops inside the
main loop:

E = 180.0
angle = numpy.min(90.0 -
numpy.arctan((win[R+1:2*R][R]-win[R][R])/dist[R+1:2*R][R]))
if angle<E: E=angle

Though maybe it will be slower for small R. You need to profile it.

Another thing not related to vectorization is that python almost
doesn't do optimizations you would expect if this code was written in
C++ or Fortran:
1. Common subexpression elimination:
win[R][R] is evaluated many times in your code, assign it to temp
variable before usage.
R + 1, R - 1, 2 * R should be computed before all the loops.

2. to get the value of numpy.arctan python has to do two dictionary
lookup in your hottest loops. Assign it to local variable at the top of
your function:
from numpy import arctan
If you are going to use min from my example, import it as well.

As regards array.array, it seems a bit complicated to reduce all my 2D
arrays to 1D arrays because I am working with the gdal library, which
works only with 'numeric' arrays.

You can try to convert numeric array into array.array.

Serge.