Getting average elevation in floating point format around NumPy point coordinate?

Question

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?

• PolyGeo: As per the 2-minute Tour, which is designed to introduce all users to this site and its protocols, there should be only one question asked per question.
• Fezter: Please edit your question to contain only one question.
• Zoltan: A similar question can be found here: gis.stackexchange.com/questions/17432/…

Answer 1

The ReadAsArray method takes starting point (xoff and yoff) and window size (xsize and ysize) arguments. To get a window around your point, you need to shift your starting point up and left and increase your size by the amount of pixels you want to buffer by.

To get a 3x3 window centred on point px, py instead of ReadAsArray(px,py,1,1), you would use ReadAsArray(px-1,py-1,3,3).

If your DEM is in floating point format, ReadAsArray will return a floating point numpy array. Alternatively, you can use numpy to cast to a different datatype e.g. vals.astype(numpy.float32) or pass in a buffer object that is an existing float numpy array.

Here's some code that will get the min, max and mean value for a window around each point (note it's a little more complicated than the example above as it handles points that fall outside the raster and windows that only partially overlap):

from osgeo import gdal
import numpy as np
if __name__ == "__main__":
 src_filename = r"C:\Temp\dem_wgs84.tif"
 #nodes = numpy.genfromtxt('tc_outlet_cor.txt')
 ## For testing
 nodes = np.array([( 148.22607422, -35.53820801),
 ( 148.25512695, -35.51940918),
 ( 148.28967285, -35.54510498),
 ( 148.26470947, -35.54309082),
 ( 154.26470947, -15.54309082)])
 # The "buffer" window in pixels, not map coordinates.
 # Note this is an N x N window, not X|Y +- N
 # For example, shape = (3, 3) is a 3x3 window which equals X|Y - 1, X|Y, X|Y + 1
 # not X|Y - 3, X|Y, X|Y + 3
 shape = (3, 3)
 offsets = [(s // 2) for s in shape]
 # Alternatively you could specify the buffer (offsets) directly instead
 # offsets = (1, 1)
 # shape = [(o*2+1) for o in offsets]
 src_ds=gdal.Open(src_filename)
 rb = src_ds.GetRasterBand(1)
 gt=src_ds.GetGeoTransform()
 rb=src_ds.GetRasterBand(1)
 for x,y in nodes:
 # Convert from map to pixel coordinates.
 # Only works for geotransforms with no rotation.
 px = int((x - gt[0]) / gt[1]) # x pixel
 py = int((y - gt[3]) / gt[5]) # y pixel
 # Skip any points outside the raster
 if 0 > px < src_ds.RasterXSize or 0 > py < src_ds.RasterYSize:
 continue
 # Reduce windows that would be partially outside the raster extent
 xoff, yoff = max([0, px - offsets[0]]), max([0, py - offsets[1]])
 xsize, ysize = min([shape[0], src_ds.RasterXSize - xoff]), min([shape[1], src_ds.RasterYSize - yoff])
 vals = rb.ReadAsArray(xoff, yoff, xsize, ysize)
 print (vals.min(), vals.max(), vals.mean())

• Heinz: so what format ReadAsArray returns depends on what format DEM is?