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
Taras
  • 32,823
  • 4
  • 66
  • 137
Ange
  • 63
  • 3

1 Answers1

7

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.

Taras
  • 32,823
  • 4
  • 66
  • 137
J. Monticolo
  • 15,695
  • 1
  • 29
  • 64
  • 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 Dec 03 '19 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. – J. Monticolo Dec 03 '19 at 22:13
  • 1
    Not a problem. It was an excellent answer even before you added that information. Thank you. – csk Dec 04 '19 at 15:34
  • 1
    Hey, please make a Pull Request instead of hiding it away here. :) – inc42 Dec 06 '19 at 07: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. – J. Monticolo Dec 06 '19 at 09:00
  • @J.Monticolo: following python script code works perfectly in Arcgis, if you convert this script in qgis python then problem will be solved. – Kapil Dev Adhikari Jan 08 '20 at 11:25
  • @KapilDevAdhikari : I've found this one but I think behind this is Arcpy and no pure Python, but the form is inspiring. Which Python script do you refer ? – J. Monticolo Jan 08 '20 at 12:17