Splitting many polygons into (almost) equal parts using a simple script in PyQGIS, the Python Console of QGIS

Fellow researchers and open-source GIS enthusiasts,

Welcome to my blog!

I’d like to start with a disclaimer – I may be a researcher of this very area but that doesn’t mean everything I do or write here will work for you, in your own desktop configurations and package versions. I have no responsibility if you lose data or mess up your installation. I also do not authorize any copies of my content, except for the gists, which are available under the MIT License.

Whatever is your geo-thing, you may have come across a Polygon which needed to be split in half, three parts, five parts (a hundred parts?). Possible applications of the split polygons are performing zonal statistics and terrain parceling (agriculture).

But wait! There is already a known pipeline used for splitting polygons in GIS Software. It was developed by Darafei Praliaskouski for PostGIS and tested by Paul Ramsey (Thank you!), and is available at Paul Ramsey’s website. The same pipeline was also adapted by Ujaval Gandhi to be used in QGIS (through the User Interface) and can be found in his blog. If you are dealing with a small number of polygons to be divided in QGIS, the aforementioned post from Ujaval Gandhi may be useful to you, especially if you are not familiar with the Python Console of QGIS (PyQGIS).

If you are a Portuguese speaker, you can also check this video on the YouTube channel of Geoaplicada, in which Jocilene Barros uses the Polygon Divider plugin.

By the way, if you are a Portuguese speaker, did you know you could read this exact blog post, and all others, in Portuguese? Just change the language by clicking the globe symbol on the upper right.

Why PyQGIS?

Oh, so many options! If you are still here, I guess you want to use PyQGIS for this task. But why would you want that? I listed four motives below:

1. You have many Polygons to divide in equal parts each. Doing the same process one-by-one is such a waste of time. Honestly, this was my motive.

2. You don’t want to install more plugins, but you also don’t want to spend too much time splitting polygons.

3. You are a Python novice/intermediate/advanced and want to see a working example of how to use PyQGIS.

4. You want to use only native tools of QGIS. This pipeline here uses only native tools.

A Disclaimer about the usage of the word “Equal”

The resulting polygons generated by pipelines which involve K-means clustering and Voronoi/Thiessen polygons (such as mine) most of the times do not generate a perfect division in terms of area. There is an error associated to each division, which we will be exploring later, on the section “Checking the results and associated errors”. Do not use my pipeline if you need perfectly equal resulting polygons.

My pipeline

For this development, I use a slightly different pipeline from the one described in Paul Ramsey’s website, but I was inspired by that pipeline to build this one.

My pipeline, in general lines:

1. Preprocessing: Split the original vector file into N other files, each one with one polygon.

2. In PyQGIS: Get the extents of the vector layer

3. Generate regular spacing points within those extents

4. Clip the regular points to the polygon area

5. Cluster your points by using K-Means algorithm, use the number of clusters as equal to the number of divisions

6. Calculate mean coordinates for each cluster of points

7. Calculate Voronoi polygons for the mean coordinates

8. Clip the polygon based on the Voronoi polygons

9. Postprocessing: check the results

The example

I start with this Shapefile:

This is its attribute table:

And, in this example, I want to divide each of the polygons in three equal parts.

I am using QGIS 3.22 Bialowieza for this tutorial.

Preprocessing details

To make sure the polygons are divided in equal (or almost equal) parts, it is important to know the area of each of them beforehand. On the attribute table, open the field calculator and create a new field called area. I gave my field the type “Decimal number (real)”. The expression is simply “$area”. Generate this column, save, and turn off the editing mode. My attribute table now looks like this:

Important: use a projected coordinate system when performing area calculations.

Another important preprocessing step. To run the PyQGIS code on this post, I need to have each polygon in a different shapefile. If you have all the polygons you want to split in the same shapefile, use “Split Vector Layer” (native:splitvectorlayer) to separate the polygons. In my case, I used the id column to split them.

Save them in a permanent folder, preferentially.

The resulting shapefiles are:

Running the PyQGIS code

This is the gist, hosted on GitHub. Made available under MIT License.