1

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?

  • 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. – PolyGeo Oct 31 '16 at 5:32
  • 1
    Please edit your question to contain only one question. – Fezter Oct 31 '16 at 5:32
  • A similar question can be found here: gis.stackexchange.com/questions/17432/… – Zoltan Oct 31 '16 at 7:58
2

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())
  • so what format ReadAsArray returns depends on what format DEM is? – Heinz Nov 1 '16 at 5:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.