2

I have a raster (geotiff) and I would like simply to know some parcentile values, specifically the 80th, 85th, 90th and 95th, because I need to highlight all values above these thresholds in QGIS.

First I thought that in Propreties--> Style-->Cumulative count cut I had simply to write the percent value in Min or Max boxes (80 and 95 respectively, for example). But I discovered that if I write 80 in the first box or in the second one (next to the cumulative count cut), the value will be different.enter image description here

So I tried to convert the raster in an Ascii grid, I opened the Ascii with text editor, I copied the cell values and I pasted in an Excel file. I deleted all the no data values (-9999 in my case) and I used the Percentile function in Excel to extract from the resulting matrix my values. I thought that it would the right way to find the values, but I wanted to be sure, so I used the Max function in Excel just to see if the value was the same to the one I can find in Propreties-->Metadata-->Maximum value.... and the to values were different, so I don't trust in what I did in Excel.

I think that what I need is a very simple opration, but I can't figure how to do it. Any help? I hope it is all clear

2

You can use pyqgis to convert your raster to numpy Array then find percentiles with numpy.percentile:

import gdal
from osgeo import gdal_array
import numpy as np

rasterfile = r"C:\Test\nh_66_5.tif" #Change

#Or replace line above with this, with the raster added to the map and activated by clicking on layer name:
#layer = iface.activeLayer()
#rasterfile = layer.source()

percentiles = [50,80,90] #Change

rasterArray = gdal_array.LoadFile(rasterfile) #Read raster as numpy array

for p in percentiles:
    print('{0}th percentile is: {1}'.format(p, np.percentile(rasterArray,p)))

enter image description here

If your raster is multiband, np.percentile will flatten the Array (=include all bands). If you want percentiles for specific band, use indexing, for example np.percentile(rasterArray[0],50)

  • Very interesting. I don't know how to use python and if I try to use the script this is what happens: rasterfile = r"Macintosh HD/Utenti/Lorenzo/Scrivania/MAPPA_FINALE_senza_Na.tif" percentiles = [50,80,90] rasterArray = gdal_array.LoadFile(rasterfile) Traceback (most recent call last): File "<input>", line 1, in <module> File "/Library/Frameworks/GDAL.framework/Versions/2.1/Python/2.7/site-packages/osgeo/gdal_array.py", line 211, in LoadFile ds = gdal.Open( filename ) RuntimeError: Macintosh HD/Utenti/Lorenzo/Scrivania/MAPPA_FINALE_senza_Na.tif: No such file or directory – Franza Mar 8 at 11:25
  • @Franza something is incorrect with the path to your raster. I dont know anything about Mac paths.. I have added some alternative lines which uses layer to get path, you could try this. – BERA Mar 8 at 11:30

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.