I am not 100% sure but it seems like gdal_rasterize takes value of first found feature by default. Is there a way to tell I would like to get minimum / maximum / average / mean (/ count / sum just to get really fancy) of vector features values within cell?

Particularly I have a csv/shapefile point file with X,Y,Z columns and points were often two or more points fall into one pixel. How can I rasterize maximum value of those within (or minimum, average, mean)?

There are sql "-where" and "-sql select_statement" but not sure if/how these could be used for this case.

  • Not with that function, see gdal.org/gdal__alg_8h.html#a50caf4bc34703f0bcf515ecbe5061a0a about the class invoked. There is no option for first/last/min/max etc.. you would need to implement your own. For points it's not that hard but you'd have to pick a language to be more specific. – Michael Stimson Aug 8 '16 at 2:18
  • @MichaelMiles-Stimson are you sure? So my best bet is to go with SAGA rasterize then? Because that has option for selection first/last/minimum/maximum/mean... – Miro Aug 8 '16 at 3:14
  • I've not used that utility so I can't say for sure. The GDAL rasterize program (which all others are based on) doesn't allow options so it would be impossible to use that utility with options without rewriting it completely - when creating an API you can always supply default values for options (thus skipping them) that exist but cannot create options that don't exist. – Michael Stimson Aug 8 '16 at 3:45

You can use GDAL Python to rasterize, however, before it is necessary a little bit of Python (in my case PyQGIS). I tried out my approach creating a little raster of 3x3 (to facilitate results validation) and, afterward, I used 'Create grid' processing tool for getting a Grid shapefile perfectly aligned with raster layer. Another shapefile, with a variable quantity of points per raster cells, was also created. This point shapefile has a field (named 'field'; see 'attribute' method in below code) with random float values. All set can be observed at the next image:

enter image description here

Next code select values in 'field' field of point shapefile and, if they are within of each feature of Grid shapefile, they are grouped by row and column indexes (for working in the case of random points with no consecutive id by cell). These values are also printed at the Python Console of QGIS.

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

points = [ feat for feat in layers[0].getFeatures() ]

grid = [ feat for feat in layers[1].getFeatures() ]

rextent = layers[2].extent()
rprovider = layers[2].dataProvider()

xmin = rextent.xMinimum()
ymax = rextent.yMaximum()
xsize = layers[2].rasterUnitsPerPixelX()
ysize = layers[2].rasterUnitsPerPixelY()

elements = []

for i, point in enumerate(points):
    posX = point.geometry().asPoint()[0]
    posY = point.geometry().asPoint()[1]

    row = int((ymax - posY)/ysize)
    col = int((posX - xmin)/xsize)

    for j, item in enumerate(grid):
        if point.geometry().within(item.geometry()):
            elements.append([row, col, point.attribute('field')])

elements.sort( key=lambda x: (x[0], x[1]) )

n = len(elements)

values = [ [] for i in range(len(grid)) ]

k = 0


for i in range(n-1):
    if elements[i][0] == elements[i+1][0] and elements[i][1] == elements[i+1][1]:
        values[k].append(elements[i + 1][2])
        k += 1

print values

After running the code at the Python Console of QGIS the printed values were:

[[58.1, 62.66, 12.93], [40.23], [93.21, 41.74], [3.23], [63.45, 56.86, 63.22, 57.62], [22.0], [48.23, 95.04], [50.23, 35.13], [43.48]]

There are nine lists because raster is 3 x 3. Each internal list has point values per each cell and its quantity is equal to number points. Now, it's easy to extract minimum, maximum or mean per cell and rasterize these values by using GDAL python module (after reshape and convert them in array with numpy).

  • Nice one, thank you. I think in my particular case it is probably more efficient to use SAGA tool Rasterize which has minimum, maximum and mean and can be called through subprocess from python the same way as gdal_rasterize. But this code gives way more flexibility, thanks again. – Miro Aug 10 '16 at 23:06

Similar problem where I had zero to many points within each grid-cell and I want the sum of an attribute for all features. I used gdal_rasterize at the CLI, (gdal v2.1.2). Standard options result in a band of all zeros.

Key was to set the initialization constant to zero and then pass the -add flag so that values from each feature were accumulated.

gdal_rasterize -te <xmin ymin xmax ymax> -tr <xres yres> -of Float32 -ot GTiff -a <targ_attrib> -init 0.0 -add -l <src_layer> src_layer.shp out_grid.tif

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.