0

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? enter image description here

  • 2
    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
  • Please read gis.stackexchange.com/help/someone-answers – user2856 Jan 6 '18 at 8:42
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.

enter image description here

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

1

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

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.