How to find buffer radius around raster cell that covers a defined sum of raster cells?

Question

I want to do a neigbourhood analysis on some raster data. I need to know how big a circle around each raster cell must be to include a certain amount of specified raster cells.

Let's say that for both rasters I want to know the radius of the circle needed to reach an aggregated value of 5. In the first case the "circle" must have a radius of 2 to reach the sum of the desired 5 because a radius of 1 would only include one raster cell that I am interested in. In the second case already a smaller circle with a radius of 1 is sufficient. My outcome should thus be a raster file where every cell has the value of the radius that is needed to reach that defined sum.

I was thinking a long time about a good title for this question. I would welcome better suggestions.

Answer 1

The answer to your question depends on the scale of your raster and if you want a unique raster "buffer" for the whole layer or per grid cell. For the first choice: The best way on QGIS might be to use the r.neighbors GRASS plugin inside the Processing toolbox for QGIS 2.0.
See a the manual here: http://grass.osgeo.org/grass65/manuals/r.neighbors.html
To do what you want you might have to build a little model in Processing, but maybe there is an easier way. The general Idea:

1. Insert your raster, start the r.neighbors tool and with a size of 3 and obviously sum as method.
2. The generated resulting raster will contain the sum of the neighborhood for each grid cell in a 3x3 distance matrix. Use the QGIS Rastercalculator to extract all cell smaller than 5 (like this RasterLayerName@1 < 5).
3. Count the name of cells inside the resulting raster (using raster stats). There is a function somewhere in Processing as well.
4. If the number of cells is equal to zero you succeeded and have found your global adequate buffer size. If not than go back to Step 1 and increase the size

Answer 2

You can also do it using grass' r.mfilter (I think).

The following filter will create a binary raster (of ones and zeros). It will produce 'one' if there are at least 5 cells in range of one cell around the center cell otherwise it will produce 0

MATRIX    3 
1 1 1 
1 0 1
1 1 1
DIVISOR   5
TYPE      P

You can re-run the function again, with the following filter, but this time it will check for cells in range of two cells.

MATRIX    5 
1 1 1 1 1 
1 0 0 0 1
1 0 0 0 1
1 1 1 1 1
DIVISOR   5
TYPE      P

Answer 3

Using GDAL/Python is another way to do this procedure. For testing it, I prepared a raster (with aleatory values of 0 and 1) and 29 rows by 29 columns. This was used for running the following code in the Python Console of QGIS (and this raster as active layer):

from osgeo import gdal
import struct
import numpy as np
layer = iface.activeLayer()
provider = layer.dataProvider()
fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 'Int32':'i', 'Float32':'f', 'Float64':'d'}
path= provider.dataSourceUri()
dataset = gdal.Open(path)
band = dataset.GetRasterBand(1)

BandType = gdal.GetDataTypeName(band.DataType)

array = []

for i in range(band.YSize - 2):
    for j in range(band.XSize - 2):
        scanline = band.ReadRaster(i, j, 3, 3, 3, 3, band.DataType)
        values = struct.unpack(fmttypes[BandType] * 9, scanline)
        tmp = np.sum(values) - values[4]
        array.append(tmp)

dataset = None

The algorithm starts at the position (x_off,y_off) = (0,0) and explores all cells considering blocks of 3x3. The element value[4] = value[2][2] is always eliminated of each sum. The neighborhood analysis is stored in the list array and could be used for creating a new raster. The complete serie obtained here (array) for this raster is exposed to continuation and it was verified that the expected values were produced.

[6.0, 4.0, 3.0, 5.0, 4.0, 6.0, 3.0, 3.0, 4.0, 2.0, 3.0, 2.0, 3.0, 4.0, 6.0, 5.0, 5.0, 3.0, 3.0, 2.0, 3.0, 4.0, 2.0, 4.0, 1.0, 3.0, 3.0, 4.0, 4.0, 5.0, 5.0, 3.0, 5.0, 4.0, 3.0, 4.0, 2.0, 4.0, 1.0, 4.0, 3.0, 4.0, 4.0, 5.0, 4.0, 5.0, 4.0, 3.0, 4.0, 4.0, 4.0, 3.0, 3.0, 3.0, 4.0, 4.0, 4.0, 5.0, 6.0, 5.0, 4.0, 4.0, 4.0, 3.0, 4.0, 3.0, 4.0, 4.0, 5.0, 5.0, 6.0, 5.0, 4.0, 2.0, 5.0, 4.0, 5.0, 4.0, 4.0, 2.0, 3.0, 3.0, 4.0, 6.0, 6.0, 5.0, 5.0, 5.0, 4.0, 3.0, 5.0, 4.0, 4.0, 3.0, 3.0, 5.0, 6.0, 5.0, 3.0, 3.0, 1.0, 4.0, 3.0, 6.0, 5.0, 4.0, 3.0, 4.0, 4.0, 4.0, 5.0, 6.0, 5.0, 6.0, 5.0, 7.0, 5.0, 4.0, 4.0, 4.0, 6.0, 5.0, 6.0, 6.0, 6.0, 6.0, 5.0, 4.0, 4.0, 4.0, 6.0, 6.0, 6.0, 5.0, 3.0, 3.0, 5.0, 3.0, 4.0, 1.0, 4.0, 3.0, 5.0, 3.0, 4.0, 4.0, 2.0, 5.0, 5.0, 8.0, 8.0, 6.0, 4.0, 3.0, 4.0, 3.0, 4.0, 4.0, 6.0, 5.0, 5.0, 5.0, 2.0, 4.0, 2.0, 3.0, 2.0, 5.0, 5.0, 7.0, 4.0, 4.0, 4.0, 4.0, 5.0, 4.0, 5.0, 5.0, 5.0, 6.0, 7.0, 7.0, 6.0, 6.0, 4.0, 6.0, 6.0, 7.0, 4.0, 5.0, 4.0, 1.0, 2.0, 3.0, 5.0, 6.0, 6.0, 6.0, 4.0, 4.0, 3.0, 4.0, 4.0, 5.0, 4.0, 4.0, 3.0, 5.0, 7.0, 7.0, 6.0, 5.0, 5.0, 3.0, 4.0, 5.0, 4.0, 4.0, 5.0, 5.0, 4.0, 5.0, 6.0, 5.0, 5.0, 4.0, 7.0, 5.0, 5.0, 2.0, 3.0, 2.0, 3.0, 5.0, 6.0, 7.0, 7.0, 6.0, 5.0, 5.0, 5.0, 5.0, 4.0, 6.0, 5.0, 5.0, 5.0, 4.0, 5.0, 4.0, 4.0, 5.0, 4.0, 3.0, 2.0, 3.0, 2.0, 2.0, 4.0, 3.0, 5.0, 4.0, 4.0, 5.0, 5.0, 6.0, 6.0, 5.0, 4.0, 6.0, 4.0, 5.0, 6.0, 7.0, 3.0, 3.0, 3.0, 5.0, 4.0, 5.0, 5.0, 4.0, 3.0, 2.0, 2.0, 3.0, 2.0, 4.0, 4.0, 3.0, 3.0, 4.0, 6.0, 7.0, 8.0, 6.0, 7.0, 3.0, 5.0, 4.0, 4.0, 4.0, 4.0, 2.0, 3.0, 4.0, 4.0, 3.0, 1.0, 2.0, 0.0, 1.0, 1.0, 2.0, 3.0, 2.0, 3.0, 3.0, 3.0, 5.0, 6.0, 7.0, 6.0, 6.0, 3.0, 4.0, 6.0, 5.0, 3.0, 2.0, 4.0, 4.0, 6.0, 5.0, 7.0, 5.0, 4.0, 3.0, 2.0, 1.0, 0.0, 1.0, 2.0, 2.0, 2.0, 3.0, 2.0, 5.0, 5.0, 6.0, 6.0, 3.0, 2.0, 3.0, 3.0, 5.0, 5.0, 4.0, 2.0, 5.0, 4.0, 5.0, 3.0, 3.0, 3.0, 1.0, 1.0, 0.0, 1.0, 2.0, 4.0, 3.0, 3.0, 2.0, 5.0, 4.0, 5.0, 3.0, 3.0, 4.0, 6.0, 6.0, 4.0, 4.0, 3.0, 5.0, 4.0, 5.0, 7.0, 6.0, 6.0, 3.0, 3.0, 2.0, 1.0, 1.0, 1.0, 4.0, 2.0, 2.0, 1.0, 3.0, 3.0, 4.0, 3.0, 5.0, 4.0, 5.0, 5.0, 5.0, 6.0, 3.0, 4.0, 2.0, 4.0, 4.0, 4.0, 4.0, 3.0, 4.0, 2.0, 2.0, 1.0, 3.0, 5.0, 5.0, 5.0, 3.0, 4.0, 3.0, 3.0, 3.0, 3.0, 5.0, 7.0, 6.0, 5.0, 3.0, 4.0, 5.0, 3.0, 3.0, 3.0, 4.0, 4.0, 3.0, 3.0, 3.0, 3.0, 1.0, 4.0, 5.0, 6.0, 5.0, 4.0, 5.0, 5.0, 3.0, 3.0, 5.0, 5.0, 4.0, 3.0, 3.0, 2.0, 2.0, 2.0, 1.0, 1.0, 1.0, 2.0, 1.0, 3.0, 4.0, 4.0, 3.0, 3.0, 4.0, 5.0, 7.0, 6.0, 7.0, 6.0, 5.0, 3.0, 3.0, 3.0, 5.0, 4.0, 3.0, 1.0, 1.0, 3.0, 5.0, 4.0, 3.0, 1.0, 3.0, 2.0, 4.0, 3.0, 6.0, 3.0, 4.0, 4.0, 4.0, 6.0, 4.0, 7.0, 7.0, 7.0, 6.0, 4.0, 5.0, 5.0, 3.0, 3.0, 2.0, 3.0, 4.0, 3.0, 2.0, 3.0, 1.0, 4.0, 2.0, 6.0, 4.0, 5.0, 2.0, 3.0, 3.0, 3.0, 5.0, 3.0, 5.0, 5.0, 6.0, 5.0, 3.0, 6.0, 5.0, 4.0, 5.0, 3.0, 3.0, 4.0, 5.0, 5.0, 4.0, 2.0, 5.0, 4.0, 6.0, 5.0, 6.0, 3.0, 3.0, 3.0, 5.0, 5.0, 5.0, 4.0, 5.0, 6.0, 7.0, 5.0, 5.0, 4.0, 4.0, 3.0, 4.0, 6.0, 7.0, 6.0, 4.0, 4.0, 4.0, 5.0, 7.0, 6.0, 4.0, 4.0, 2.0, 4.0, 3.0, 4.0, 4.0, 3.0, 1.0, 4.0, 4.0, 5.0, 2.0, 6.0, 3.0, 3.0, 3.0, 3.0, 2.0, 3.0, 5.0, 4.0, 5.0, 3.0, 4.0, 5.0, 4.0, 5.0, 4.0, 4.0, 4.0, 6.0, 4.0, 6.0, 4.0, 4.0, 5.0, 4.0, 6.0, 4.0, 3.0, 4.0, 3.0, 4.0, 3.0, 4.0, 4.0, 5.0, 3.0, 5.0, 3.0, 5.0, 6.0, 4.0, 4.0, 2.0, 6.0, 4.0, 5.0, 1.0, 4.0, 2.0, 4.0, 4.0, 3.0, 5.0, 3.0, 2.0, 3.0, 3.0, 4.0, 3.0, 4.0, 3.0, 5.0, 5.0, 5.0, 4.0, 3.0, 2.0, 2.0, 4.0, 3.0, 5.0, 4.0, 6.0, 4.0, 6.0, 4.0, 4.0, 5.0, 4.0, 6.0, 5.0, 2.0, 3.0, 5.0, 4.0, 6.0, 4.0, 5.0, 5.0, 5.0, 4.0, 2.0, 3.0, 4.0, 5.0, 5.0, 4.0, 6.0, 5.0, 4.0, 3.0, 3.0, 4.0, 5.0, 6.0, 5.0, 6.0, 6.0, 2.0, 1.0, 3.0, 1.0, 4.0, 3.0, 6.0, 6.0, 5.0, 6.0, 3.0, 4.0, 3.0, 4.0, 4.0, 3.0, 3.0, 4.0, 4.0, 6.0, 6.0, 5.0, 3.0, 4.0, 6.0, 6.0, 5.0]

A new evaluation, with a slight modification of the above code, could be executed with 5x5 blocks if the requeriment of sum >= 5 is not accomplished for 3x3 blocks .