From the polygonsplitter plugin, I adapted the essential part of the code for QGIS 3.x :
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) ] )
#QgsMapLayerRegistry.instance().addMapLayer(polyLayer)
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)
#Figure out direction to start with from cutting method
#If alternating, start horizontally
if method=="a":
firstDirection="h"
else:
firstDirection=method
self.alternatingSlice(geom,polyLayer,targetArea,granularity,firstDirection,method)
else:
self.addGeomToLayer(geom,polyLayer)
polyLayer.updateExtents()
#if self.debug: print recurs
QgsProject.instance().addMapLayer(polyLayer)
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())))
self.addGeomToLayer(multiGeom,polyLayer)
#self.addGeomToLayer(QgsGeometry.fromPolyline([startPt,endPt]),lineLayer)
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:
self.addGeomToLayer(geom,polyLayer)
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
def addGeomToLayer(self,geom,layer):
"""
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().addFeatures([fet])
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, select your layer to cut in the layer tree under QGIS, and execute the following code :
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.