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?

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>'

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.