# How to split polygon into equal area polygons in QGIS 3?

There is a plugin in QGIS 2.18 called `'Polygon Divider'` which divides an irregular polygon into smaller equal-area polygons of a specified/custom-sized but in QGIS 3.x version there is no `'Polygon Divider'` plugin or tool to perform that task.

So, how to divide/split the irregular polygon into equal area/parts polygons in QGIS 3.6? Is there any Python script to split such type of polygon?

As shown in the below picture:

• I don't know this plugin. But, you can try to create a grid first (processing - toolbox - create Grid - grid type : [] Polygon) , and than clip this layer with your polygon's layers. – Christophe P. Apr 27 '19 at 7:56
• You could upgrade the 2.x plugin to work in 3.x github.com/qgis/QGIS/wiki/Plugin-migration-to-QGIS-3 – Ian Turton Apr 27 '19 at 8:19
• @Cyril I have updated my question and added screenshot too. – Kapil Dev Adhikari Dec 22 '19 at 16:17
• I added my answer. The study was made using an experimental script written in the programming language pgSQL/PostGIS SQL... – Cyril Mikhalchenko May 12 '20 at 16:02
• Another important detail: "Do not believe your eyes..." – Cyril Mikhalchenko May 29 '20 at 8:17

All of the necessary tools can be found from processing toolbox and you have to be a bit innovative to achieve this. Basically by adopting this method you could achieve pretty good results, although it doesn't give you 90 degree angles in all places, but roughly equal area:

1.Create enough random points inside a polygon. I did 10 000

1. Run k-means clustering for the random points. Use as many clusters as many areas you want:

1. Run Vector Geometry --> Aggregate and use the cluster id as the aggregation field. This way you will get a multipoint layer from the clusters.

2. Run Vector geometry --> centroids for the aggregated points

3. Run Voronoi polygons for the centroids. Use enough buffer to cover the whole original polygon. For me this was around 40%

4. Create lines from the voronoi polygons (Vector geometry --> polygons to lines). Currently I think you can only split a polygon layer with lines and not with another polygon layer (not 100% sure about this, but this works for sure).

5. Finally run Split with Lines tool and use your original polygon layer as input and voronoi lines as the split layer:

In this attempt, I tried to accomplish @Christophe's idea suggested in comments. Please don't be harsh, it is definitely askew and fragile nevertheless it is somehow working. The script does not include the debugging part.

Proceed with `Plugins > Python Console > Show Editor` and paste the script below

``````# defining required imports
import sys

# defining a class
class LayerIntoParts:

# initiating layer's name to split and number of parts provided from user
def __init__(self, layer_name, segments):
self.layer_name = layer_name
self.segments = segments

# setting up the input layer and checking it's type
def layer_geometry_check(self):
layer = QgsProject.instance().mapLayersByName(self.layer_name)[0]
if not layer.isValid():
sys.exit

# checking geometry type, will work only with polygons and multipolygons
if layer.wkbType() == QgsWkbTypes.Polygon:
layer_valid = layer
# converting multipart into single part
elif layer.wkbType() == QgsWkbTypes.MultiPolygon:
layer_valid = processing.run('native:multiparttosingleparts',
{'INPUT': layer,
'OUTPUT': "memory:"}
)["OUTPUT"]
# giving an error for not polygons and multipolygons
else:
raise ValueError("This geometry type is not supported!\Only Polygons/Multipolygons are allowed.")
sys.exit

# returning the validated layer
return layer_valid

# defining required parameters
def layer_and_parameters(self):
# inheriting the output from 'layer_geometry_check' method
layer = self.layer_geometry_check()

# getting the layer's area
area_total = [f.geometry().area() for f in layer.getFeatures()]
area_total = area_total[0]

# getting the layer's extent and adjusting it for processing functions
layer_extent = layer.extent()
layer_extent_string = "{},{},{},{}".format(layer_extent.xMinimum(), layer_extent.xMaximum(), layer_extent.yMinimum(), layer_extent.yMaximum())

# building a dictionary to store the output parameters
param = {
'layer': layer,
'layer_crs': layer.crs(),
'layer_extent': layer_extent_string,
'layer_area_total': area_total,
'area_segment': area_total / self.segments
}

# returning the dictionary with output parameters
return param

# creating grid for the input layer
def layer_2_grid(self):
# inheriting the output from 'layer_and_parameters' method
param = self.layer_and_parameters()

# creating grid for the layer with 'qgis:creategrid' algorithm
layer_grid = processing.run('qgis:creategrid',
{'TYPE': 2,
'CRS': param['layer_crs'],
'EXTENT': param['layer_extent'],
'HSPACING': 10000,
'VSPACING': 10000,
'HOVERLAY': 0,
'VOVERLAY': 0,
'OUTPUT': "memory:"}
)["OUTPUT"]

# deleting redundant attributes (appear by default when using the algorithm to create a grid)
layer_grid_clean = processing.run('qgis:deletecolumn',
{'INPUT': layer_grid,
'COLUMN': ['left', 'right', 'top', 'bottom'],
'OUTPUT': "memory:"}
)["OUTPUT"]

# returning layer as a grid
return layer_grid_clean

# clipping the grid with the input layer
def clip_grid_and_layer(self):
# inheriting the output from 'layer_2_grid' method
grid = self.layer_2_grid()
# inheriting the output from 'layer_and_parameters' method
layer = self.layer_and_parameters()['layer']

# clipping the grid with the layer by means of 'qgis:clip' algorithm
layer_clip = processing.run('qgis:clip',
{'INPUT': grid,
'OVERLAY': layer,
'OUTPUT': "memory:"}
)["OUTPUT"]

# returning the clipped grid
return layer_clip

# categorizing the clipped grid layer
def categorize(self):
# inheriting the output from 'clip_grid_and_layer' method
layer_clip = self.clip_grid_and_layer()
# inheriting the output from 'layer_and_parameters' method
param = self.layer_and_parameters()

# transforming clipped grid layer's elements into a list with
# id, calculated area, and geometry
elements = []
for f in layer_clip.getFeatures():
case = {"f_id": f.id(), "area": f.geometry().area(), "geom": f.geometry()}
elements.append(case)

# lists for used cells' ids
used_ids_i = []
used_ids_j = []
# variable for categories
n = 0

# looping through all cells
for i in elements:
area_calc = i['area']
used_ids_i.append(i['f_id'])

for j in elements:

# conditions: calculated area less than a "golden" area of a segment, no duplicates of cells, and sufficient intersection of edges between cells
if area_calc <= param['area_segment'] and j['f_id'] not in used_ids_i and i['f_id'] not in used_ids_j and (i['geom'].intersection(i['geom'])).length() > 0:
used_ids_j.append(j['f_id'])
area_calc += j['area']
n += 1

# signing a category to each cell
i['cat'] = 'cat_' + str(n)
i['area_calc'] = area_calc

# returning categorized cells as list
return elements

# dissolving categorized cells and outputting the result
def dissolve_result(self):
# inheriting the output from 'layer_and_parameters' method
param = self.layer_and_parameters()
# inheriting the output from 'categorize' method
elements_raw = self.categorize()

# creating an empty later for output
vl = QgsVectorLayer("Polygon?crs={}".format(param['layer_crs'].authid()), "output", "memory")

# accessing the empty layer
pr = vl.dataProvider()

# embedding new field values
QgsField("cat", QVariant.String),
QgsField("area_calc", QVariant.Double),
QgsField("area_gold", QVariant.Double)
])

# updating the empty layer
vl.updateFields()

# looping through a list of categorized cells and inserting them into the empty layer
for elem in elements_raw:
feat = QgsFeature()
feat.setGeometry(elem['geom'])
feat.setAttributes([elem['f_id'], elem['cat'], round(elem['area_calc'], 8), round(param['area_segment'], 8)])

vl.updateExtents()

# dissolving by category by means of 'qgis:dissolve' algorithm
layer_dissolved = processing.run('qgis:dissolve',
{'FIELD': ['cat'],
'INPUT': vl,
'OUTPUT': 'memory:'}
)["OUTPUT"]

# deleting 'cat' attribute using 'qgis:deletecolumn' algorithm
layer_dissolved_clean = processing.run('qgis:deletecolumn',
{'INPUT': layer_dissolved,
'COLUMN': ['cat'],
'OUTPUT': "memory:"}
)["OUTPUT"]

# updating 'id' attribute
with edit(layer_dissolved_clean):
for feature in layer_dissolved_clean.getFeatures():
feature.setAttribute(feature.fieldNameIndex('id'), feature.id())
layer_dissolved_clean.updateFeature(feature)

# returning function output

my_class = LayerIntoParts('test_area_4', 2)
my_class.dissolve_result()
``````

Adjust the `my_class = LayerIntoParts('name of your layer', num_of_parts)`

Press `Run script` and get the output

Testing. I am using a polygon of Corsica that was taken from Natural Earth | Admin 0 – Countries

Splitting into 2 parts with

``````my_class = LayerIntoParts('test_area_4', 2)
my_class.dissolve_result()
``````

Splitting into 3 parts with

``````my_class = LayerIntoParts('test_area_4', 3)
my_class.dissolve_result()
``````

Notes: (TODO)

• Parameters `'HSPACING'` and `'VSPACING'` have to be defined in advance
• Does not work with Multipolygons

References:

I recently "partially" updated the Split Polygon plugin to QGIS 3, created by Uroc Ilic. It is now available in the QGIS repository.

But it works only for convex polygons, that is, not very irregular.

This is an approach that could be tried:

1) So far, it is most likely impossible to accurately separate a complex polygon by horizontal and vertical lines.

2) The problem lies in the combination of straight and curved lines.

3) It is possible to separate a complex polygon accurately only with curved lines.

4) Yes, with small errors it is possible to divide the polygon by horizontal and vertical lines and they will accumulate at the second (third, etc.) level, i.e. the 1st level will allow me to precisely divide the curvilinear polygon into 2 exactly equal parts, and from the second level, when we have sides with straight and curved lines, the geoinstrument will start to accumulate errors...

"...until the road is overwhelmed" :-).

5) An example of the study is shown in the figures below

and

6) Maybe one day I'll reconsider my point of view, but that's what came to me now :-)...

Next step to think about is how to find the average value of the polygon and then bring the area of sub-plots to the average value of the polygon.

Now, there is a ‘Polygon Divider’ plugin for QGIS 3.x version, which is the solution.

• Please edit your question to include a detailed description of the steps necessary to accomplish the tasks. Otherwise, it is likely that the link will be obsolete in the future and, as a result, the answer becomes obsolete. – Aaron Sep 27 '20 at 20:42