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

Question

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?

Answer 1

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.

Answer 2

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

def circularmask(r=5):

    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)
    array[mask] = True

    return array

if __name__ == '__main__':

    radius = 5
    mask = circularmask(radius)
    xsize, ysize = mask.shape

    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

        xoff, yoff = (px - radius, py + radius)

        array = rb.ReadAsArray(px, py, xsize, ysize)
        print(zonalmean(array, mask, nodata))