5

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?

map for the treatment block


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
5
  • 3
    Are they regular rectangles ? Commented Dec 3, 2019 at 15:04
  • No they are irregular
    – Ange
    Commented Dec 3, 2019 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
    Commented Dec 3, 2019 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%. Commented Dec 3, 2019 at 15:53
  • 2
    gis.stackexchange.com/a/333756/120129 Commented Dec 3, 2019 at 15:55

1 Answer 1

7
+50

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) ] )
            #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, copy the code above in a tab in the QGIS Python Code Editor, execute it with the Run script 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 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.

7
  • 2
    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
    Commented Dec 3, 2019 at 18:26
  • 1
    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. Commented Dec 3, 2019 at 22:13
  • 1
    Not a problem. It was an excellent answer even before you added that information. Thank you.
    – csk
    Commented Dec 4, 2019 at 15:34
  • 1
    Hey, please make a Pull Request instead of hiding it away here. :)
    – inc42
    Commented Dec 6, 2019 at 7:30
  • 1
    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. Commented Dec 6, 2019 at 9:00

Not the answer you're looking for? Browse other questions tagged or ask your own question.