I am working on a script that automatically generates heatmaps from point data. I am aware of the heatmap/KDE tool already existing in PyQGIS (I am working in QGIS 3.10). What I want to do after computing the heatmap is to reclassify the data into five different classes like so:

enter image description here

Instead of this:

enter image description here

To do so, a reclassify is required. This is doable with the reclassify tool in QGIS, however, you must insert the reclass values yourself. This is not desirable when creating a 'one size fits all' script, as one would need to generate the min and max values and intervals for each class over and over again. In ArcPy, the slice tool can 'slice' the heatmap into different sections according to for example the 'natural breaks' standard.

Is there a good way to do so in PyQGIS?

Note: The points where the heatmap is created with include values, which should be taken into account when creating the heatmap.

ps. I got it to work with quantiles with the help of numpy. The result looks like this:

import numpy as np 
from osgeo import gdal, gdal_array

dataset = gdal.Open("path/to/rasterlayer.sdat")
array = dataset.ReadAsArray()

output = np.zeros_like(array).astype(np.uint8)

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)

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 = ("path/to/rasterlayer_percentiles.tif")
gdal_array.SaveArray(output, outname, "gtiff", prototype=dataset)

Would there be a similar way to do this, but then with natural breaks?

1 Answer 1


To do Jenks breaks directly, import a suitable Python library. Here's an example of some overlapping normal distributions.

from jenkspy import jenks_breaks
import numpy as np
import matplotlib.pyplot as plt

d = np.concatenate([a,b,c])

histdata = plt.hist(d,100)
breaks = jenks_breaks(d,3)

Results in :

[-2.3007428299285735, 2.944479747738755, 8.264521464618356, 15.028677383261819]

which seems correct.

enter image description here

  • Wingnut - I am indeed aware of this function. The problem is however that I am automating a script for 'any' dataset with 'any' values. Therefore, a new imported table has to be written 'on the go', every time a new dataset is entered.
    – Bram
    Commented Apr 26, 2021 at 14:19
  • I see. Then you will need to load the dataset using osgeo.gdal, get the raster as Numpy, reclassify the Numpy array, and resave as a new dataset :-)
    – wingnut
    Commented Apr 26, 2021 at 14:35
  • That's what I did. I have provided the code that I used in the question itself. However, this only works for quantiles/percentiles. I was wondering if there was a way to reclassify with natural breaks/jenks within numpy or in any other way. An answer to that would be the very answer to my question! :)
    – Bram
    Commented Apr 26, 2021 at 14:39
  • I can't see an algorithm anywhere (even online), just the methodology. I would rephrase the original question and just ask if anyone knows how to compute Natural Breaks, without without Jenks optimisation. I was confused by the images, as they can be done in QGIS by changing the layer styling (without explicit computation). Hopefully, that makes sense.
    – wingnut
    Commented Apr 26, 2021 at 14:57
  • 1
    I did find this: pip install jenks - from medium.com/analytics-vidhya/… - and anaconda has jenkspy - in Jupyter, you can type jenkspy? for usage - I learned something today. Thanks.
    – wingnut
    Commented Apr 26, 2021 at 15:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.