# Natural Breaks Results Difference in Different GIS Analyst Tools

I calculated natural break clustering in each GIS software and Python Jenkspy library.

The data is population of each city in Europe and this is the result below.

Jenkspy library and QGIS shows somewhat similar results (but not really) while other trials return totally different results.

I am going to try with CartoDB PostGIS Extension though, I am already doubting the result.

the input data is all same and I made 6 classes for division. But I am not understanding why the results are different.

Does anyone know the reasons?

Surprisingly, QGIS 2.18 and QGIS 3 gave me different results with same data.

• Can you share the data or can you reproduce this with `pop_est` for all features in `ne_110m_admin_0_countries.shp`? Commented May 3, 2019 at 10:17
• I added another trials Commented May 4, 2019 at 15:04

The first question to investigate is, what method are these programs using to calculate Natural Breaks? Are they all using the same method?

So it seems that all four programs use Jenks natural breaks optimization. What's not yet answered (but you can probably find this out if you dig further into the documentation for each program) is if they use exactly the same method, or if they each use their own method based on Jenk's optimization. For example, maybe Geoda uses a combination of Jenks optimization and Fisher's discriminant analysis.

The Jenks optimization method... is a data clustering method designed to determine the best arrangement of values into different classes. (source)

Jenks natural breaks optimization can be done in one of two ways.

1. The first way is an iterative process. It creates class breaks (possibly arbitrarily, possibly by some other method), calculates the "sum of the squared deviations from the class means," checks the result against "a minimal value," and if it doesn't meet that "minimal value," it repeats the steps and checks again. This goes on until the minimal value is met.
2. The second way is to calculate all the possible combinations of class breaks, calculate "sum of the squared deviations from the class means" for all combinations, and pick the combination with the lowest value.

To summarize:

• All four programs use Jenks optimization, but there's a lot of variation in how this method can be applied, leading to different results.
• You can test which program has the "best" method. Calculate the "sum of the squared deviations from the class means" for the class breaks created by each program. Whichever program has the lowest value, is getting the closest to optimal results.
• thanks for the detail explanation! Commented May 7, 2019 at 0:49

I'm replying to this post very late so it may not be very useful, but Jenkspy (of which I'm the developer) and QGIS implement the same logic, judging by their code. Which is also the same logic implemented in MapClassify Python package as "FisherJenks"

Whether for Jenkspy, MapClassify or QGIS (and even the classInt library in R), the core of the algorithm implemented is sometime referred to as Fisher's "exact optimization" method, a reference implementation of which can be found in Fortran in cmlib: https://lib.stat.cmu.edu/cmlib/src/cluster/fish.f (see also comments in QGIS code about this).

The main difference between Jenkspy and QGIS is that QGIS takes a sample of the input dataset (cf. QGIS code: the variable `mMaximumSize` is set to 3000 in the header file, so if there is more than 3000 values in the dataset, the code samples 3000 values or 10% of the dataset, whichever is larger) whereas Jenkspy, and MapClassify (when using "FisherJenks" and not "FisherJenksSampled"), use all the values given to it as input.