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

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

# 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?

• 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
• 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

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) / gt)  # x pixel
py = int((y - gt) / gt)  # 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]), max([0, py - offsets])
xsize, ysize = min([shape, src_ds.RasterXSize - xoff]), min([shape, 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