2

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?
  • As per the Tour there should be only one question asked per question. – PolyGeo Jan 23 '18 at 8:56
2

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)
  • Right now the statistics are added to a new column in the polygon shapefile I use to query/limit the raster with (and all 558 files median value cant fit since the shape only allows 255 columns). So this is not optimal since I now have to cut the data out of the attribute table, rearrange it in excel and then join it with the second shapefile. This shapefile is a point layer containing 558 features, all named identically with the 558 rasters mentioned above. Is there a way to get this done directly in Python? Query median by polygon and then join by filename to the point shape instead? – johlund Jan 23 '18 at 16:30
  • @johlund - Most welcome, glad it helped! In regards to your question, I think that should be doable but I would suggest asking a new question and include all the relevant details and code you have tried :) – Joseph Jan 24 '18 at 10:28
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.

  • Thanks a lot xunilk! Your suggested code runs without errors but I can't find the median results anywhere? It's not appended into the polygon shape attribute that have loaded in QGIS. I don't have any problem using numpy instead of zonalstatistics, but just so that I understand correct: where in your code will my polygon shapefile narrow the area that it calculates median on? I'm asking because this numpy approach takes much longer time then zonal stats takes to calculate four values (average, rejects, max and min), so I'm just wondering if it calculates median for the complete raster? – johlund Jan 23 '18 at 16:19
  • If the above is sorted out, I would love to hear you're ideas on how to append the median values into a second shapefile (a point layer with features named just like the rasters) – johlund Jan 23 '18 at 16:20
  • What does your point vector layer represent? Centroids of each raster layer extent? – xunilk Jan 23 '18 at 17:37
  • The rasters are viewshed results from different locations, the point layer contains all the observerpoints. 558 rasters with the extent of the whole city of Stockholm and the point vector layer contains the 558 observer points. The third file, the polygon vector, is a planned new development area (much smaller extent then the rasters), and it's based on that polygon I wanna calculate median to be written to the point layer for visualisation purpouse. Rasters are named 1.tiff, 2.tiff, 3.tiff etc and the points have a column called name with values 1, 2, 3, etc. – johlund Jan 23 '18 at 17:59

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.