By Tim K

2017-02-03 13:46:08 8 Comments

The task I want to complete is essentially exactly what the Zonal Statistics plugin for QGIS aims to achieve: gives statistics about the values of a raster layer in areas defined by polygons.

However, the way Zonal Statistics seems to work is that it treats a pixel as either 'in' or 'out' of the polygon ('in' if more than half the pixel is in?). Due to limitations of my data, my raster pixel size (30m) is close to the size of the polygons I am using (100m radius circles) so those edge effects become very important. I would therefore like an alternative method which calculates the same statistics, but based on the actual area of each pixel contained within the polygon.

One option I have thought of would be to crop the raster to match my polygons, polygonize the raster, export the values and areas of the generated polygons, and carry out the calculations separately, but that seems like a lot of effort compared to a tool like Zonal Statistics. This method may also be complicated by the fact that some of my polygons overlap.


@xunilk 2017-02-06 21:19:47

Zonal Statistics, which include partial pixels, can be performed by using PyQGIS. I tried out PyQGIS code with a small raster (20 rows x 20 columns, pixel size 30 m) and a shapefile with two features that are corresponding to 100 m radius circles; as it can also be observed at next image. However, only the first feature was used in the calculations.

enter image description here

Next code determines the four vertex points for each cell; but it's only enough that one of them is within each circle to get complete cell geometry. Cell geometry is intercepted with circles features and its centroid is used to get raster values by using 'identify' method of QgsRasterDataProvider class. The indx_inter is the overlapping index based in areas of intersected zone and cells. It can be used as weighting factor for each recovered pixel value.

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

feat_buffer = layers[1].getFeatures().next() #first feature of buffer layer

prov_raster = layers[2].dataProvider()
raster_extent = layers[2].extent()
xmin, ymin, xmax, ymax = raster_extent.toRectF().getCoords()

rows = layers[2].height()
cols = layers[2].width()

sizeX = layers[2].rasterUnitsPerPixelX()
sizeY = layers[2].rasterUnitsPerPixelY()

points = []

x = xmin

for i in range(rows):
    for j in range(cols):
        points.append([QgsPoint(x + i*sizeX, ymax + j*(-sizeY)),
                       QgsPoint(x + (i+1)*sizeX, ymax + j*(-sizeY)),
                       QgsPoint(x + (i+1)*sizeX, ymax + (j+1)*(-sizeY)),
                       QgsPoint(x + i*sizeX, ymax + (j+1)*(-sizeY))])
    x = xmin

for item in points:
    geom1 = feat_buffer.geometry()
    if QgsGeometry.fromPoint(item[0]).within(geom1) or \
       QgsGeometry.fromPoint(item[1]).within(geom1) or \
       QgsGeometry.fromPoint(item[2]).within(geom1) or \
           geom2 = QgsGeometry.fromPolygon([item]) #cell
           geom3 = geom2.intersection(geom1) #intersection cell and buffer
           indx_inter = geom3.area()/geom2.area() 
           centroid = geom3.centroid().asPoint()
           value = prov_raster.identify(centroid,
           print indx_inter, value

After running the code, it's printed simultaneously at the Python Console of QGIS indx_inter and value (pixel value) values. They can be combined to compute weighting zonal statistics (where indx_inter equal 1 is corresponding to an intersected area identical to cell geometry).

0.0625204194487 5.0
0.233286907572 5.0
0.128285659439 7.0
0.430554583203 9.0
0.964620626089 7.0
1.0 3.0
0.998752162581 8.0
0.670804885299 6.0
0.0382110109824 6.0
0.259108201581 10.0
0.997982751942 8.0
1.0 2.0
1.0 7.0
1.0 7.0
1.0 6.0
0.569683803093 7.0
0.610928539309 8.0
1.0 2.0
1.0 10.0
1.0 9.0
1.0 2.0
1.0 7.0
0.923521388879 8.0
0.636179042117 9.0
1.0 4.0
1.0 2.0
1.0 2.0
1.0 10.0
1.0 1.0
0.948771891687 10.0
0.341485379809 4.0
1.0 3.0
1.0 8.0
1.0 4.0
1.0 4.0
1.0 8.0
0.654078229379 7.0
0.00328203975005 6.0
0.592744891916 1.0
0.999302262394 3.0
1.0 6.0
1.0 3.0
0.817116132863 4.0
0.0908059107671 3.0
0.192046339914 3.0
0.397494464342 5.0
0.29124537879 6.0
0.0178963092058 9.0

Last polygon (indx_inter = 0.0178963092058; partial intersection) and centroid of iteration were visualized by using QuickWKT plugin. Raster value was recovered with Value Tool plugin (see next image).

enter image description here

@rumski20 2017-02-08 14:32:27

I like the novel approach to solving the problem, especially the complete code example and accompanying visuals. One thing to note is that this approach should not be used with rasters containing categorical data (ie. land cover).

@rumski20 2017-02-03 14:10:01

I think you hit on your best option, which is to convert the raster to a vector and then intersect the result with your polygon layer.

As a way of explanation regarding the frustration you're experiencing with Zonal Statistics (and actually this will also apply to your idea of cropping or "clipping" the raster), there is no alternative way for this to operate; a single raster cell has to be either in or out.

You also mention that some of your "zones" (ie. polygons) overlap. This will also negate results of a standard Zonal Statistics tool since most of those explicitly mention that they do not account for this and recommend running the tool over each feature separately if this is the case.

Related Questions

Sponsored Content

1 Answered Questions

[SOLVED] Zonal stats for complete pixels qgis

2 Answered Questions

0 Answered Questions

Zonal Statistics calculator differences between QGIS Lyon and Wien

0 Answered Questions

QGIS Zonal Statistics for polygons is missing results?

0 Answered Questions

Zonal statistics behavior

1 Answered Questions

1 Answered Questions

1 Answered Questions

[SOLVED] Zonal Statistics Freezing QGIS 2.2

2 Answered Questions

[SOLVED] How to Calculate Custom Zonal Statistics?

Sponsored Content