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[0]
res = []
i = 0
while i < rows:
mx,my=points[i][0], points[i][1] #coord in map units
#Convert from map to pixel coordinates.
#Only works for geotransforms with no rotation.
px = int((mx - gt[0]) / gt[1]) #x pixel
py = int((my - gt[3]) / gt[5]) #y pixel
intval=rb.ReadAsArray(px,py,1,1,buf_type=gdal.GDT_Float64)
# print intval[0] #intval is a tuple, length=1 as we only asked for 1 pixel value
res.append(intval[0])
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?