There are a couple Stackexchange posts asking this question, but neither of them were sufficiently answered:

I am writing Python scripts that perform essentially the same operations as Zonal Stats, and I have been using Zonal Stats to check my scripts' outputs.

I have noticed that when there are many cells in the polygon being queried, my results match perfectly, but there are discrepancies when there are only a few raster cells in the polygon. The following is an example of this. I am computing "count" for each polygon.

Exhibit 1:

enter image description here The returned count is 0.45, which looks roughly correct if pixels can be treated fractionally rather than simply in/out.

Exhibit 2:

enter image description here Similar result; count is 0.39.

Exhibit 3:

enter image description here

With a larger poylgon, now we see the count is exactly 2. Note that there are exactly 2 pixel centers contained in the polygon.

Exhibit 4:

enter image description here

Here, only one pixel center is within the polygon. The count is 1.008, though. If you exclude the polygon area from the pixel whose center is in the polygon, the remaining area is obviously greater than 0.008 of a pixel.

From these exhibits, it seems that Zonal Statistics treats these cases specially. It looks like to me that if there are less than 2 pixel centers within the polygon, it performs some kind of averaging such as polygon_area/raster_cell_area to return the count. I should note that I have 48 of these polygons, and all returned counts greater than or equal to 2 are integers; fractional counts are only returned for <2.

How are these fractional-pixel counts incorporated into the other statistics output by Zonal Stats?

How does Zonal Statistics work exactly? indicates strange behavior with max/min when using large raster cells, and my own testing returns mean values not congruous with simple "pixel center in/out of the polygon." To reiterate, I get exact agreement when the underlying raster contains many pixels within the polygon. I couldn't find any documentation about this, and it may be important when people have never had problems with Zonal Stats but use a coarse raster and remain unaware that there is a different behavior when there are only a few pixels within the polygon.

up vote 1 down vote accepted

QGIS first makes an initial pass, checking to see if the center of each raster cell is within the polygon. If fewer than two cell centers are within the polygon, it performs a vector-based intersection for all intersecting cells, whether their center is within the polygon or not, and computes a weight that is the fraction of each cell that is covered by the polygon. If two or more cell centers are within the polygon, those cells are assigned a weight of 1 and all other cells are assigned a weight of 0.

How these weights are applied depends on the statistic in question.

  • The sum is calculated as the sum of each pixel value times its weight.
  • The count is calculated as the sum of all pixel weights.
  • The mean is calculated as the sum divided by the count.
  • The median, variance, standard deviation, minimum, and maximum are calculated for all pixels with a weight > 0, but do not take individual pixel weights into account in the calculations.

Relevant source code:

Selection of pixel inclusion/exclusion algorithm

Statistics calculation

I can handle partially overlapping pixels adapting code (for PyQGIS in QGIS 3) in the answer to this question: Robust Zonal stats in R and QGIS.

import processing

registry = QgsProject.instance()

Polygons = registry.mapLayersByName('Polygons_new')

iterator = Polygons[0].getFeatures()

feats_Polygons = [ feat for feat in iterator ]

Raster = registry.mapLayersByName('Raster')
rprovider = Raster[0].dataProvider()

xmin, ymin, xmax, ymax = Raster[0].extent().toRectF().getCoords()

extent = str(xmin)+ ',' + str(xmax)+ ',' +str(ymin)+ ',' +str(ymax) 

parameters = {"TYPE":2, 
              "EXTENT":extent, 
              "HSPACING":1000, 
              "VSPACING":1000,
              "HOVERLAY":0,
              "VOVERLAY":0, 
              "CRS":Raster[0].crs().authid(), 
              "OUTPUT":"/home/zeito/pyqgis_data/tmp_grid.shp"}

path = processing.run('qgis:creategrid', 
                         parameters)

grid = QgsVectorLayer(path['OUTPUT'],
                      'grid',
                      'ogr')

feats_grid = [ feat for feat in grid.getFeatures() ]

geom_pol = []
attr = []

for featg in feats_grid:
    for featp in feats_Polygons:
        if featg.geometry().intersects(featp.geometry()):
            geom = featg.geometry().intersection(featp.geometry())
            geom_pol.append(geom.asWkt())
            pt = geom.centroid().asPoint()
            value = rprovider.identify(pt,
                                       QgsRaster.IdentifyFormatValue).results()[1]
            attr.append([value, geom.area()])

epsg = Polygons[0].crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer&field=value:double&field=area:double&field=w_value:double""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'polygon',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(geom_pol)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i, attr[i][0], attr[i][1],(attr[i][0]*(attr[i][1]/1e6))])
    feat.setGeometry(QgsGeometry.fromWkt(geom_pol[i]))

prov.addFeatures(feats)

QgsProject.instance().addMapLayer(mem_layer)

Polygons_new referred at above code can be observed at following image:

enter image description here

After running it at Python Console of QGIS 3, result was:

enter image description here

Attributes table has weighted values for pixel for each intersected polygon area (with a temporal grid that coincides with raster cells).

At following example can be observed a shapefile (Polygons_new2) that intersects a complete pixel. Its weighted value at attributes table is the same as pixel value (as expected). Remaining pixel values are proportional to covered area by shapefile.

enter image description here

  • Thanks for the contribution; I see you make a polygon grid from your raster grid to find the relative contribution of each pixel. This would be nice to incorporate into Zonal Statistics to handle these coarse raster cases. The question about Zonal Statistics's behavior remains. – Jon Mar 23 at 14:57

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.