# How to get mean values of adjacent cells at point location when extracting values to points?

I use this code extract raster values to points.

``````from osgeo import gdal,ogr
import struct

src_filename = '/tmp/test.tif'
shp_filename = '/tmp/test.shp'

src_ds=gdal.Open(src_filename)
gt=src_ds.GetGeoTransform()
rb=src_ds.GetRasterBand(1)

ds=ogr.Open(shp_filename)
lyr=ds.GetLayer()
for feat in lyr:
geom = feat.GetGeometryRef()
mx,my=geom.GetX(), geom.GetY()  #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

structval=rb.ReadRaster(px,py,1,1,buf_type=gdal.GDT_UInt16) #Assumes 16 bit int aka 'short'
intval = struct.unpack('h' , structval) #use the 'short' format code (2 bytes) not int (4 bytes)

print intval[0] #intval is a tuple, length=1 as we only asked for 1 pixel value
``````

Base on this, I want get mean cells values adjacent at the points. Such as showing in the graph: the value of point is the mean of valid pixels values in the red circle. I assume the radius is n pixels. Nodata values should be ignored. Are there any existing function for this processing?

• One way would be to buffer the Points and then use Zonal Statistics – BERA Jul 11 '17 at 11:58
• or compute the mean for each pixel in raster and extract it – radouxju Jul 11 '17 at 12:25
• As pointed out @BERA, you only need to create a buffer point layer for your points and to add a few lines in your code to do that. – xunilk Jul 11 '17 at 14:12
• Have you tried anything at all? That code is unchanged from the other answer. – user2856 Jul 11 '17 at 23:53

As pointed out @BERA, you only need a few lines in your code to do that. Code was modified as follow (paths are for my example):

``````from osgeo import gdal,ogr
import struct
from qgis.analysis import QgsZonalStatistics

src_filename = '/home/zeito/pyqgis_data/utah_demUTM2.tif'
shp_filename = '/home/zeito/pyqgis_data/random_points.shp'

layer = iface.activeLayer()  #this is the buffer point layer

zoneStat = QgsZonalStatistics(layer,
src_filename,
"",
1,
QgsZonalStatistics.Mean)

zoneStat.calculateStatistics(None)

src_ds=gdal.Open(src_filename)
gt=src_ds.GetGeoTransform()
rb=src_ds.GetRasterBand(1)

ds=ogr.Open(shp_filename)
lyr=ds.GetLayer()
for feat in lyr:
geom = feat.GetGeometryRef()
mx,my=geom.GetX(), geom.GetY()  #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

structval=rb.ReadRaster(px,py,1,1,buf_type=gdal.GDT_UInt16) #Assumes 16 bit int aka 'short'
intval = struct.unpack('h' , structval) #use the 'short' format code (2 bytes) not int (4 bytes)

print intval[0] #intval is a tuple, length=1 as we only asked for 1 pixel value
``````

After creating a buffer point layer (and did it as active layer) for my points, I ran above code and got a field named 'mean' in its attributes table; as expected.

To consider radius as n pixels, to be used as distance buffer, you need first to calculate diagonal of cell raster.

Here's another way that doesn't use QGIS. It simply builds a circular mask using `numpy.ogrid` (modified from this SO answer), then calculates the mean from the masked data. Note this does not handle edge cases (where radius > distance to the edge of the raster).

``````from osgeo import gdal,ogr
import numpy as np

def nodata_to_nan(array, nodata):
array = array.astype(np.float64)
array[array == nodata] = np.nan
return array

def zonalmean(array, zone, nodata=None):
if nodata is not None:
return np.nanmean(nodata_to_nan(array, nodata)[zone])
else:
return(array[zone].mean())

n = r*2+1
y,x = np.ogrid[-r:r+1, -r:r+1]
mask = x*x + y*y <= r*r

array = np.zeros((n, n), dtype=np.bool)

return array

if __name__ == '__main__':

src_filename = 'test.tif'
shp_filename = 'test.shp'

src_ds=gdal.Open(src_filename)
gt=src_ds.GetGeoTransform()
rb=src_ds.GetRasterBand(1)
nodata = rb.GetNoDataValue()

ds=ogr.Open(shp_filename)
lyr=ds.GetLayer()
for feat in lyr:
geom = feat.GetGeometryRef()
mx,my=geom.GetX(), geom.GetY()  #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