# Split polygon in equal parts in QGIS [duplicate]

I have a polygon of an area of `46,465m²` and a perimeter of `136,1m`.

It is a rectangular polygon and I want to split it into three equal parts (which are my treatment blocks).

Could anybody give me some ideas on how to proceed, please?

I am trying to run the code, however, it gives me some error when I introduce my target area, which it is `15.487m²` (the total area of my polygon is `46,465m²`). Why is that?

``````Slicing, No of part:  1

Slicing, Granularity remaining:  30

Traceback (most recent call last):
File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\code.py", line 90, in runcode
exec(code, self.locals)

File "<input>", line 1, in <module>

File "<string>", line 6, in <module>

File "<string>", line 45, in splitSelected

File "<string>", line 86, in alternatingSlice

File "<string>", line 136, in cutPoly

File "<string>", line 204, in signedDistCentroidFromLine

NameError: name 'math' is not defined
``````
• Are they regular rectangles ? – J. Monticolo Dec 3 '19 at 15:04
• No they are irregular – Ange Dec 3 '19 at 15:09
• The plugins that they suggest to install do not exist in the QGIS 3.4.I am not able to solve the problem by following those suggestions – Ange Dec 3 '19 at 15:45
• Your problem is recurrent and difficult to automate. If you only have to do this task once, the easiest way is to do it by hand. You create two lines inside and then split. You calculate the area of the 3 polygons and adjust as you go if you want to get a percentage of 33.3%. – Vincent Bré Dec 3 '19 at 15:53

From the polygonsplitter plugin, I adapted the essential part of the code for QGIS 3.x :

``````import math

class EqSplitPolygon:
#def __init__(self,iface):
def __init__(self):
self.debug=True
pass

def splitSelected(self,targetArea,granulFactor,method="h",splitEven=True):
global recurs
recurs=0;
layer = iface.mapCanvas().currentLayer()
if layer:
#Gets layer CRS for new layer
crs=layer.crs().description()
if self.debug: print("Starting, Layer crs: " + crs)
# Create a new memory layer and add an area attribute
polyLayer = QgsVectorLayer("MultiPolygon?crs="+crs, "split_poly", "memory")
polyLayer.dataProvider().addAttributes( [ QgsField("area", QVariant.Double) ] )
allFeatures=False
if not layer.selectedFeatures():
layer.invertSelection();
allFeatures=True
#save original target area
origTargetArea=targetArea
# Loop though all the selected features
for feature in layer.selectedFeatures():
geom = feature.geometry()
if self.debug: print("Starting Number of original geoms: ", str(len(geom.asGeometryCollection())))
if self.debug: print("Starting Number of part to split into: ", str(geom.area()/targetArea))
div=round(geom.area()/origTargetArea)
if div<1:
div=1
if splitEven:
targetArea=geom.area()/div
if self.debug: print("Spliteven selected. modifying target area to:", targetArea)
if div>1:
granularity=round(granulFactor*geom.area()/targetArea)
if self.debug: print("Granularity: ", granularity)
#If alternating, start horizontally
if method=="a":
firstDirection="h"
else:
firstDirection=method
self.alternatingSlice(geom,polyLayer,targetArea,granularity,firstDirection,method)
else:
polyLayer.updateExtents()
#if self.debug: print recurs
if allFeatures:
layer.invertSelection();

def alternatingSlice(self,geom,polyLayer,targetArea,granularity,direction,method):
"""
Slice a poly in alternating directions
"""
global recurs
recurs+=1
if self.debug: print("******************************")
if self.debug: print("Slicing, No of part: ", str(recurs))
if self.debug: print("Slicing, Granularity remaining: ", str(granularity))
bbox=[
geom.boundingBox().xMinimum(),
geom.boundingBox().yMinimum(),
geom.boundingBox().xMaximum(),
geom.boundingBox().yMaximum()
]
if direction=="h":
step=(bbox[2]-bbox[0])/granularity
pointer=bbox[0]
else:
step=(bbox[3]-bbox[1])/granularity
pointer=bbox[1]
totalArea=0
slices=0
#save the original geom
tempGeom=QgsGeometry(geom)
#start slicing until targetArea is reached
while totalArea<targetArea*0.999:
pointer+=step
if direction=="h":
startPt=QgsPointXY(pointer,bbox[1])
endPt=QgsPointXY(pointer,bbox[3])
(multiGeom,tempGeom)=self.cutPoly(tempGeom,startPt,endPt)
else:
startPt=QgsPointXY(bbox[0],pointer)
endPt=QgsPointXY(bbox[2],pointer)
(tempGeom,multiGeom)=self.cutPoly(tempGeom,startPt,endPt)
if multiGeom!=None:
totalArea+=multiGeom.area();
slices+=1
if self.debug: print("Slicing, Slices: ", str(slices))
#do the real cutting when reached targetArea and add "left" feature to layer
if self.debug: print("Cutting with line, Cutline:", startPt,",",endPt)
if direction=="h":
(multiGeom,geom)=self.cutPoly(geom,startPt,endPt,True)
if multiGeom:
if self.debug: print("After split, Parts to the left:", str(len(multiGeom.asGeometryCollection())))
if geom:
if self.debug: print("After split, Parts to the right:", str(len(geom.asGeometryCollection())))
else:
(geom,multiGeom)=self.cutPoly(geom,startPt,endPt,True)
if geom:
if self.debug: print("After split, Parts above:", str(len(geom.asGeometryCollection())))
if multiGeom:
if self.debug: print("After split, Parts under:", str(len(multiGeom.asGeometryCollection())))
if geom:
if geom.area()>targetArea:
if (method=="v") or ((method=="a") and (direction=="h")):
self.alternatingSlice(geom,polyLayer,targetArea,granularity-slices,"v",method)
else:
self.alternatingSlice(geom,polyLayer,targetArea,granularity-slices,"h",method)
else:

def cutPoly(self,geom,startPt,endPt,debug=False):
"""
Cut a geometry by a 2 point line
return geoms left of line and right of line
"""
#if we have disjoint Multi geometry as geom to split we need to iterate over its parts
splittedGeoms=[]
leftFragments=[]
rightFragments=[]
#if self.debug: print "Number of geoms when slicing: ",str(len(geom.asGeometryCollection()))
for geomPart in geom.asGeometryCollection():
#split the actual part by cut line defined by startPt,endPt
(res,splittedGeomsPart,topo)=geomPart.splitGeometry([startPt,endPt],False)
splittedGeoms+=splittedGeomsPart
#Add the remaining geomPart to the rightFragments or letfFragments
#depending on distance
d=self.signedDistCentroidFromLine(geomPart,startPt,endPt)
if d>0:
rightFragments.append(geomPart)
else:
leftFragments.append(geomPart)
#if self.debug: print j,splittedGeoms

for fragment in splittedGeoms:
"""
calculate signed distance of centroid of fragment and the splitline
if signed distance is below zero, the point is to the left of the line
if above zero the point is to the right of the line
"""
d=self.signedDistCentroidFromLine(fragment,startPt,endPt)
#if debug==True:
#if self.debug: print d

if d > 0:
rightFragments.append(fragment)
else:
leftFragments.append(fragment)

#if self.debug: print "Left frags:",len(leftFragments),"Right frags:",len(rightFragments)
leftGeom=self.buildMultiPolygon(leftFragments)
rightGeom=self.buildMultiPolygon(rightFragments)
return leftGeom,rightGeom

def buildMultiPolygon(self,polygonList):
"""
Build multi polygon feature from a list of polygons
"""
geomlist=[]
for geom in polygonList:
# Cut 'MULTIPOLYGON(*) if we got one'
if geom.asWkt()[:12]=="MULTIPOLYGON":
geomWkt=geom.asWkt()[13:len(geom.asWkt())-1]
else:
# Cut 'POLYGON' if we got one
geomWkt=geom.asWkt()[7:]
geomlist.append(str(geomWkt))
multiGeomWKT="MULTIPOLYGON("
multiGeomWKT +=",".join(geomlist)
multiGeomWKT+=")"
#if self.debug: print multiGeomWKT
multiGeom=QgsGeometry.fromWkt(multiGeomWKT)
return multiGeom

"""
Add a geometry to a layer as a new feature
No attributes are set
"""
fet = QgsFeature()
fet.setGeometry(geom)
area=geom.area()#/1000000
if self.debug: print("Area of geom added to layer:", str(area))
layer.dataProvider().changeAttributeValues({fet.id(): { 0: area}});
layer.updateExtents()

def signedDistCentroidFromLine(self,geom,startPt,endPt):
#calculate signed distance of centroid of fragment and the splitline
v1=endPt[0]-startPt[0]
v2=endPt[1]-startPt[1]
A=v2
B=-v1
C=-v2*startPt[0]+v1*startPt[1]
centr=geom.centroid().boundingBox()
return (A*centr.xMinimum()+B*centr.yMinimum()+C)/math.sqrt(A**2+B**2)
``````

For use this class, copy the code above in a tab in the QGIS Python Code Editor, execute it with the `Run script` . Add a tab with the green +, copy the code below, select your layer to cut in the layer tree under QGIS, modify the `targetArea`, `granulFactor` and execute the following code again with `Run script` :

``````eqsplit_inst = EqSplitPolygon()
eqsplit_inst.splitSelected(
targetArea= 123.456,  # area of the polygon / 3
granulFactor=10,  # higher the number is, more precise is the cut
method="h",  # h for vertical, v for horizontal
splitEven=True
)
``````

The result isn't exact but works for simple cuts.

• For those of us who only ever took one coding class ten years ago in a different coding language (I think it was Java, but I'm not 100% sure. There were lots of brackets.), can you explain what to do with this code? Would you paste it into the "custom function" tab of the expression builder? Or do you paste it into the python console? A screenshot with the right button/location in QGIS circled would be super helpful. – csk Dec 3 '19 at 18:26
• Sorry, I was in a rush when I post it, I just wanted to provide the QGIS 3.x version because it was just few edits. I added some precisions for use the scripts. – J. Monticolo Dec 3 '19 at 22:13
• Not a problem. It was an excellent answer even before you added that information. Thank you. – csk Dec 4 '19 at 15:34
• Hey, please make a Pull Request instead of hiding it away here. :) – inc42 Dec 6 '19 at 7:30
• a PR need to adapt the whole plugin code, test it, and it will be hidden if the maintainer don't accept it. But I'll consider it. – J. Monticolo Dec 6 '19 at 9:00