How does QGIS Zonal Statistics handle partially overlapping pixels?

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:

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

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.