7

I have a raster file which has a range (min: 80 max: 120). I would like to reclass this raster by using quantiles. Is that possible?

7

If it is just for visualisation, then you can adjust how the raster is displayed in QGIS, by choosing single band pseudo colour with discrete colour interpretation - you can then adjust the boundaries between each colour yourself.

Raster display in qgis

If you actually need to produce an output raster classified by percentile, then it may or may not be possible to do so with the raster calculator in QGIS. But I would strongly advise that you have a look at the python tools that are available, as they really are ideal for this. Numpy has a percentile function which does exactly this.

I believe that QGIS comes with numpy as standard, so you can do this in the QGIS python console (open with ctl + alt + p).

An example script would be

import numpy as np
from osgeo import gdal, gdal_array

# open the dataset and retrieve raster data as an array
dataset = gdal.Open("/path/to/image.tif")
array = dataset.ReadAsArray()

# create an array of zeros the same shape as the input array
output = np.zeros_like(array).astype(np.uint8)

# use the numpy percentile function to calculate percentile thresholds
percentile_80 = np.percentile(array, 80)
percentile_60 = np.percentile(array, 60)
percentile_40 = np.percentile(array, 40)
percentile_20 = np.percentile(array, 20)
percentile_0 = np.percentile(array, 0)

The numpy.where function to change the output array; the syntax is np.where((condition), x, y)), where x is the value to set if the condition evaluates to true at that index, and y is the value to set if the condition evaluates to false.

output = np.where((array > percentile_0), 1, output)
output = np.where((array > percentile_20), 2, output)
output = np.where((array > percentile_40), 3, output)
output = np.where((array > percentile_60), 4, output)
output = np.where((array > percentile_80), 5, output)

The gdal_array.SaveArray is a very handy function which allows you to specify a prototype, in this case the input dataset, from which the projection info is copied.

outname = "/path/to/output_name.tif")
gdal_array.SaveArray(output, outname, "gtiff", prototype=dataset)

Documentation for numpy.percentile

| improve this answer | |
  • No, it is not only for visualisation. I will try the options that I got from this forum. – Pimpel Feb 25 '17 at 19:45
  • Ah ok, I wasn't sure about that. Let me know if you have any trouble with the python method, you should be able to just copy and paste it into the python console in qgis, after changing the file names and percentile values. – user6072577 Feb 25 '17 at 20:47
4

You can do in QGIS processing toolbox. Use the model builder and build a model with r.quantile followed by r.recode. Something like this: enter image description here

| improve this answer | |
  • The model could be useful in some cases, but it will put the brake value into the cell instead of using a categorical (integer) value. – César Arquero Feb 2 '18 at 8:45
0

@Hossein Madadi thanks for writing useful code, but it would not tackle the "nan" values and do not gives quantile as it supposed to be, so I have modified the code and its working like charm :)

import numpy as np
from osgeo import gdal, gdal_array
import glob
import subprocess
import os

base_path = r'D:\etc\\'

# open the dataset and retrieve raster data as an array
input_raster = base_path + "Other\\rainfall_idw.tif"
dataset = gdal.Open(input_raster)
band = dataset.GetRasterBand(1)
array = band.ReadAsArray()

nodata_val = band.GetNoDataValue()
if nodata_val is not None:
    array = np.ma.masked_equal(array, nodata_val)

array_ignored_nan = array[array >= array.min()]

# create an array of zeros the same shape as the input array
output = np.zeros_like(array).astype(np.uint8)

# use the numpy percentile function to calculate percentile thresholds
percentile_80 = np.percentile(array_ignored_nan, 80)
percentile_60 = np.percentile(array_ignored_nan, 60)
percentile_40 = np.percentile(array_ignored_nan, 40)
percentile_20 = np.percentile(array_ignored_nan, 20)
percentile_0 = np.percentile(array_ignored_nan, 0)

print(percentile_0, percentile_20, percentile_40, percentile_60, percentile_80)

output = np.where((array > percentile_0), 1, output)
output = np.where((array > percentile_20), 2, output)
output = np.where((array > percentile_40), 3, output)
output = np.where((array > percentile_60), 4, output)
output = np.where((array > percentile_80), 5, output)

outname = os.path.splitext(input_raster)[0] + "_Classified.tif"
gdal_array.SaveArray(output, outname, "gtiff", prototype=dataset)

This will produce raster, with nodata value equal to zero, because classified raster is "uint", and it will show as black, so in order to make it transparent you need to add transparency in produced raster, following is code for that,

Shapefile = base_path + 'shps\\study_area.shp'
Shapefile = Shapefile.replace("\\\\", "\\")
Image = outname.replace("\\\\", "\\")

OutImage = (OutDir + os.path.split(Image)[1].split('.')[0] + '_Clipped.tif').replace("\\\\", "\\") # Defines Output Image

if os.path.exists(OutImage):
    os.remove(OutImage)
    
# Clip image
subprocess.call('gdalwarp -srcnodata "0" -dstnodata "9999" -cutline "' + Shapefile + '" -crop_to_cutline -dstalpha "' + Image + '" "' + OutImage + '"', shell=True, stdout=open(os.devnull, 'w'), stderr=subprocess.STDOUT)
| improve this answer | |

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.