I just started working with QGIS, and for the life of me I cannot figure out how to solve this problem:
I have several points of interest, these are stored in a shapefile (just simple points). I have (for now) one raster file with values. What I want to do basically is select those cells in my raster that are near the points in the shapefile (I want several window-sizes, like first only one cell, than 9 (3x3), than 25 (5x5) etc). I want to average these cells and get the median and stdev.
I want to do this in a Python script, since I will probably have to do this again (and again) at some point. I am new to QGIS, but in ArcGIS I would have tried converting my points to a raster, and 'grow' this raster while calculating zonal statistics.
However since in QGIS it should be easier to get data-level access, I was hoping to avoid all that and just convert my points to a position in a raster, and than loop over the values using Python. I already figured out how to get to the data using ReadAsArray:
raster = gdal.Open(inputRaster) band = raster.GetRasterBand(1) data = band.ReadAsArray()
However I am not sure how (if it is even possible once the data is converted to an array) to link the positions in my shapefile to the position in the array.
If this is not possible is there another strategy I should consider?