Storing result from gdallocationinfo as variable in Python

Question

I'm trying to find the cell value at various locations in a set of rasters so I can count the number of real values (aka not 'NoData') that occur at all the locations within each raster. The various locations are represented as points in a point shapefile. I exported the shapefile's attribute table (including lat/long) to a csv and am looping over each point to get the value. I have a script ready to go, but am having trouble storing the output of the 'gdallocationinfo' as a result.

How can I do this in python?

Here is my pseudo code:

for raster in rasters:
    count = 0
    for i in range(0, Npts):
        lat = points[i][10] # lat is 10th row of csv
        lon = points[i][9] # longitude is 9th row of csv

        result = os.system('gdallocationinfo -valonly -wgs84 %s %s %s' % (raster, lon, lat))
        print result

        if not result == 'nan': count+=1 # if result is no 'nan', add to the count

    row = [raster, str(count)]
    # write row to csv

When I run the actual code, the system call for 'gdallocationinfo' prints out the correct values (or 'nan'). But when it gets to 'print result', it just prints 0. So the command call is giving me the correct result, but this obviously is not the way to store the result as a variable so I can use it later. I'm calling the python code from the linux command line.

Answer 1

You can capture standard out this way:

from subprocess import Popen, PIPE

cmd = ['gdallocationinfo', '-geoloc', '-valonly', src, lng, lat]
p = Popen(cmd, stdout=PIPE)
p.wait()
val = p.stdout.read()

Answer 2

I did a little more digging and found the answer to my question:

instead of using os.system, the correct syntax for storing the result in a variable is:

result = os.popen('gdallocationinfo -valonly -wgs84 %s %s' % (lyr, loc)).read()

Answer 3

I think you will have to use the subprocess syntax, it is explained in this post:

subprocess to call gdal from within python

Answer 4

To speed up the process you could use gdallocationinfo with input/output redirection. Export all your lat/lon values to input.csv and read the corresponding cell values from output.csv:

os.system(r'gdallocationinfo -valonly -wgs84 "%s" <%s >%s' % ('dhm.tif','input.csv','output.csv'))
with open('output.csv','r') as f:
    result = f.read()