Created
June 17, 2020 13:40
-
-
Save danzig666/12df413e22f47ff8b2ca167973f2ba1a to your computer and use it in GitHub Desktop.
Polygon Splitter QGIS script
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| """ | |
| Splitting polygons into equal area parts | |
| converted to qgis3 from https://github.com/csandor/polygonsplitter | |
| """ | |
| #from __future__ import print_function | |
| #from builtins import str | |
| #from builtins import object | |
| from qgis.core import * #(QgsFeature, QgsGeometry, | |
| #QgsVectorLayer, QgsMapLayerRegistry, | |
| #QgsField,/) | |
| from PyQt5.QtCore import * #from qgis.PyQt.QtCore import QVariant | |
| from qgis.utils import iface | |
| import math | |
| def printl(message): | |
| if 'sself' in globals() and sself is not None: | |
| sself.dockwidget.label.append(message) | |
| else: | |
| print(message) | |
| QApplication.processEvents() | |
| #this is required for the qgis script runner plugin | |
| def run_script(iface): | |
| #global lineLayer | |
| #lineLayer = QgsVectorLayer("LineString?crs=3785", "split_line", "memory") | |
| eqsplit_inst = EqSplitPolygon(iface) | |
| eqsplit_inst.splitSelected(15,5,"v",True) | |
| #QgsMapLayerRegistry.instance().addMapLayer(lineLayer) | |
| class EqSplitPolygon(object): | |
| #def __init__(self,iface): | |
| def __init__(self): | |
| self.debug=False | |
| pass | |
| def splitSelected(self,targetArea,granulFactor,method="h",splitEven=True): | |
| global recurs | |
| recurs=0; | |
| layer = iface.layerTreeView().currentLayer() #iface.mapCanvas().currentLayer() | |
| if layer: | |
| #Gets layer CRS for new layer | |
| crs=layer.crs().description() | |
| if self.debug: # fix_print_with_import | |
| printl("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: # fix_print_with_import | |
| print("Starting Number of original geoms: ",str(len(geom.asGeometryCollection()))) | |
| if self.debug: # fix_print_with_import | |
| 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: # fix_print_with_import | |
| print("Spliteven selected. modifying target area to:", targetArea) | |
| if div>1: | |
| granularity=round(granulFactor*geom.area()/targetArea) | |
| if self.debug: # fix_print_with_import | |
| 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) #QgsMapLayerRegistry.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: # fix_print_with_import | |
| print("******************************") | |
| if self.debug: # fix_print_with_import | |
| print("Slicing, No of part: ",str(recurs)) | |
| if self.debug: # fix_print_with_import | |
| 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: # fix_print_with_import | |
| print("Slicing, Slices: ", str(slices)) | |
| #do the real cutting when reached targetArea and add "left" feature to layer | |
| if self.debug: # fix_print_with_import | |
| print("Cutting with line, Cutline:", startPt,",",endPt) | |
| if direction=="h": | |
| (multiGeom,geom)=self.cutPoly(geom,startPt,endPt,True) | |
| if multiGeom: | |
| if self.debug: # fix_print_with_import | |
| print("After split, Parts to the left:",str(len(multiGeom.asGeometryCollection()))) | |
| if geom: | |
| if self.debug: # fix_print_with_import | |
| 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: # fix_print_with_import | |
| print("After split, Parts above:",str(len(geom.asGeometryCollection()))) | |
| if multiGeom: | |
| if self.debug: # fix_print_with_import | |
| 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.exportToWkt())-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: # fix_print_with_import | |
| 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) | |
| eqsplit_inst = EqSplitPolygon() | |
| eqsplit_inst.splitSelected(40000,5,"a",True) # area 40000, granulFactor 5, method: "a": alternating, "h":horizontal, "v": vertical, splitEven: True |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment