1

What I'm trying to do:

  1. Extract median by a single polygon from multiple GeoTIFF rasters stored in a folder (Solved)
  2. Add each median value to a point shapefile (Match/join by: id column for each point in shape = raster filename)

Basically the project has three sources:

  1. One folder containing 558 GeoTIFF raster's which are viewshed results from different observer points. All raster's have the same extent and are named 1.tiff, 2.tiff, 3.tiff .... 558.tiff

  2. One vector point layer with the above observer points. The ID field is numbered from 1 to 558 (which corresponds to the raster filenames)

  3. One vector polygon layer which represents a planned new development area (much smaller extent then the raster's). This is the polygon that raster median is calculated from.

What I got so far:

By using Python console from within Qgis and Qgs.ZonalStatistics I have, with help from @Joseph, been able to extract median values by polygon (3) from all TIFF's in the folder (1) with this code:

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 new planned development polygon layer I use to query the raster's with (3). The columns are named by filename with "median" added without speces behind the .tiff extention, see example below:

enter image description here

But 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 observer points vector layer, and on top of that I have to do it in two batches since 255 is the column limit in a shapefile.

How can I re-write the above code to take the median output I get from Qgs.ZonalStatistics and append it (join by id/filename) with the vector point layer containing all the observer points (2)?

1
+50

You could try the following which when calculating the median, stores the name of the raster and the median value in a dictionary. Then we can create a single new field (e.g. "Median") in your point layer, loop through the features and if the FID matches the name of the raster, we can add the median value to a median field:

import glob, os
from qgis.analysis import QgsZonalStatistics

vectorlayer = iface.activeLayer()
rasterfolder = 'D:/effect_test/'

pointLayer = QgsMapLayerRegistry.instance().mapLayersByName("Plangrans")[0]
field_index = pointLayer.fieldNameIndex('Median')
median_dict = {}

os.chdir(rasterfolder)
for lyr in glob.glob("*.tiff"):
    raster_name = os.path.splitext(os.path.basename(lyr))[0]
    median_value = QgsZonalStatistics(vectorlayer, lyr, attributePrefix=lyr, rasterBand=1, stats=QgsZonalStatistics.Statistics(QgsZonalStatistics.Median)).calculateStatistics(None)
    median_dict[raster_name] = median_value

with edit(pointLayer):
    for feat in pointLayer.getFeatures():
        if feat['FID'] in median_dict.keys():
            pointLayer.changeAttributeValue(feat.id(), field_index, median_dict[feat['FID']])
  • 1
    You're amazing, thank you very much for helping me out, again! Haven't been able to try this yet but I'm marking it as the accepted answer based on how smooth everything went the last time you helped. Hoping that if I run into any trouble running this I could reach out and see if you might have a clue on how to solve it. I'll also get back as soon as I tried it to let you know if it worked or not. – johlund Jan 30 '18 at 8:05
  • @johlund - Most welcome, although you really should have waited to accept the answer until after you ran it! But I ran it on some examples so hopefully it will work for you too :) – Joseph Jan 30 '18 at 10:05
  • Didn't want you to think I was ungrateful! – johlund Jan 30 '18 at 10:35
  • @johlund - You said "thank you", that's more than enough for me =) – Joseph Jan 30 '18 at 10:43

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.