Based on How do I get the pixel value of a GDAL raster under an OGR point without NumPy?, I wrote the following python function to find elevation height values of given coordinates in a DEM file:
from osgeo import gdal,ogr import struct import numpy def single_point_elev(dem_file, points): src_filename = dem_file src_ds=gdal.Open(src_filename) gt=src_ds.GetGeoTransform() rb=src_ds.GetRasterBand(1) rows = points.shape res =  i = 0 while i < rows: mx,my=points[i], points[i] #coord in map units #Convert from map to pixel coordinates. #Only works for geotransforms with no rotation. px = int((mx - gt) / gt) #x pixel py = int((my - gt) / gt) #y pixel intval=rb.ReadAsArray(px,py,1,1,buf_type=gdal.GDT_Float64) # print intval #intval is a tuple, length=1 as we only asked for 1 pixel value res.append(intval) i = i + 1 nres = numpy.array(res).reshape(len(res), 1) return nres if __name__ == "__main__": nodes = numpy.genfromtxt('tc_outlet_cor.txt') print single_point_elev('hc.tif', nodes)
and the result was:
[[ 107.] [ 109.] [ 108.] [ 26.] [ 25.] [ 27.] [ 247.] [ 239.] [ 249.] [ 281.] [ 283.] [ 281.]]
In the case above, input DEM file was a region containing streams, but as I checked the result, I found that some height values which belonged to downstream section were higher than that of upstream points, and besides, I got two points had the same elevation height value, as showed in the result above.
How could I get average elevation height value around given coordinate in floating point format?