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?

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())