Create elevation-service from GeoTIFF using Python and GDAL

Question

I'd like to create a little server for us, that tells us the elevation of a point in a GeoTIFF file. I already got it working basically using the following:

dataset = gdal.Open('test_file.tif', gdal.GA_ReadOnly)
geotransform = dataset.GetGeoTransform()
px, py = convertLatLngToPixels(geotransform, lat, lng)
val = self.dataset.GetRasterBand(1).ReadAsArray(px, py, 1, 1)

Now the issues is, that our real raster file, is a VRT that combines several GeoTIFFs. Using my code now on that VRT does not work any more...

Does anybody know any solution to using a VRT file?

• user30184: Study the source code of gdal.org/programs/gdallocationinfo.html#gdallocationinfo or use that as is.
• Georg: But I'd need it for python... But i guess i could use the result of the lifonly option and then continue inside my code? But the issue will be that gdallocationinfo -lifonly is not really fast... Is there a fast way to find the file that I'd need to use inside python?
• Jon: What do you mean "does not work anymore"? Throws an error? Gives the wrong value? ... I don't see why it wouldn't work on a vrt as opposed to a tif.
• ThomasG77: Borrow code from github.com/OSGeo/gdal/blob/master/gdal/swig/python/gdal-utils/… maybe?

Answer 1

As it was said in comments, you can find the answer in gdallocationinfo sources. It does the following thing:

In [37]: band = ds.GetRasterBand(i)
In [38]: band.GetMetadataItem('Pixel_1000_1000', 'LocationInfo') 
Out[38]: '<LocationInfo><File>/files/Files/GIS/raster_nsk/test.tif</File></LocationInfo>'
In [39]: band.GetMetadataItem('Pixel_992644_2000', 'LocationInfo') 
Out[39]: '<LocationInfo><File>/files/Files/GIS/raster_nsk/abc1.tif</File><File>/files/Files/GIS/raster_nsk/abc2.tif</File></LocationInfo>'

• bugmenot123: The questioner wants pixel values, your answer seems to show something else.

Answer 2

Maybe you find an answer for your question here: Profile Tool has too high resource consumption: Alternative or work around?

In my solution I'm using a VRT file that combines ~5000 8MB TIF files to a virtual DEM mosaic for building a terrain profile, and I'm really impressed about the processing speed of Python GDAL.