Python: Extract median by polygon from multiple raster files

Question

What I'm trying to do:

• Extract median (or 50th quantile) by a single polygon from 558 GeoTIFF rasters stored in a folder

What I got so far:
Using Python from within Qgis and Qgs.ZonalStatistics I've managed to extract statistics by polygon from *.tiff files in a directory, and add each value to a new column in the above polygon shapefile (column name based on raster filename). The code looks like this:

import glob, os, qgis.analysis
vectorlayer = qgis.utils.iface.mapCanvas().currentLayer()
rasterfolder = 'D:/effect_test/'
os.chdir(rasterfolder)
for lyr in glob.glob("*.tiff"):
    qgis.analysis.QgsZonalStatistics(vectorlayer, lyr, attributePrefix=lyr, rasterBand=1).calculateStatistics(None)

What's left to resolve:

• With the above code I'm getting statistics for average, rejects, max and min while I need to extract the median (which is the same as the 50th quantile). I've been trying to add QgsZonalStatistics.Mean suggested in this question and also QgsZonalStatistics::Mean suggested in the QgsZonalStatistics documentation. But it fails telling me 'QgsZonalStatistics' is not defined. How to re-write above code to get median or the 50th quantile instead?

Answer 1

Maybe something like the following:

import glob, os
from qgis.analysis import QgsZonalStatistics
vectorlayer = iface.activeLayer()
rasterfolder = 'D:/effect_test/'
os.chdir(rasterfolder)
for lyr in glob.glob("*.tiff"):
    QgsZonalStatistics(vectorlayer, lyr, attributePrefix=lyr, rasterBand=1, stats=QgsZonalStatistics.Statistics(QgsZonalStatistics.Median)).calculateStatistics(None)

Answer 2

I read your question before being edited and reopened and it could be preferable to calculate medians for each particular raster instead to use QgsZonalStatistics. So, you can re-write your code as follows:

import numpy as np
import glob, os, qgis.analysis
rasterfolder = '/home/zeito/pyqgis_data/effect_test' #change for your folder
os.chdir(rasterfolder)

attr = []

for lyr in glob.glob("*.tif"):  #change for your tiff extension
    layer = QgsRasterLayer(lyr, lyr)
    provider = layer.dataProvider()
    extent = provider.extent()
    rows = layer.height()
    cols = layer.width()
    block = provider.block(1, extent, cols, rows)
    values = [ block.value(i,j) for i in range(rows) for j in range(cols) ]
    attr.append([layer.name(), np.median(values)])

print attr

Medians are calculated by using 'median' method from numpy.

Advantage of this approach is that later you can easily create your point layer for incorporating raster's names and calculated medians in its attributes table.