Split polygons into equal areas within each region

Question

I've a set of regions within a city. I'd like to split each region into five subregions and then get the centroid of each one of them.

I know how to make a vector grid, but the problem I have is that the generated points are not taking into account the regions, they're only related to the general map.

I have this:

And I want to split each region in five equal size areas. Something like this:

I'd like to build a visualization similar to this one: http://bl.ocks.org/mbostock/8814734. In order to accomplish that, I need to split each region, because the original ones are pretty big.

Answer 1

Ok, I had a go at this using existing tools for QGIS...

---

1. Download/enable the Buffer by percentage plugin from:

   Plugins > Manage and Install Plugins...
   

   This creates buffers for each of your features.

---

2. Run this plugin on your layer and use the Buffer area percentage option. Enter 20 (for 20%), save the output and run it again for 40, 60, and 80.

   Here is an example layer:

   

   And here is running the plugin:

   

   And here are the results of the plugin (it will help to add the percentage value in the name as we will use this to identify the layers):

   

---

3. Now we need to run this really ugly model (or you can download it and copy it to your /.qgis2/processing/models folder):

   

   When you run it, you will need to enter the correct layers according to parameter (hence why we need to include the % value in the name):

   

---

4. The output should hopefully have split your features into five equal parts. I tested this on several areas and calculated the area using the Field Calculator with the expression $area. Although it might not look attractive to what you wanted in your image, it does seem to split the features equally...

   

Answer 2

I will suggest solution in PostGIS. Use ST_Subdivide, calculate total number of vertices of each polygon and divided it by 5. Use this result as max_vertices in query. Further info can be found here. http://postgis.net/docs/manual-2.2/ST_Subdivide.html

Answer 3

I had a go at this in ArcMap since i dont know how to do it in QGIS. Maybe someone else can translate the ideas into QGIS. Im using ArcGIS 10.4 with Advanced License (Avanced is needed for Erase Tool which is used in script).

Script below will divide all polygons in an input feature class in five almost equal area parts. Almost equal area because the script does not check that all parts are exactly the same area, only that the area of each part is within a specified margin of error.

If the polygon being processed is round(ish) it will be divided like a cake:

If it's not round, it's elongated and will be divided like this:

But not like this (divisions not parallel):

Example:

Code:

import arcpy,math
arcpy.env.workspace = r'C:\TEST.gdb'
arcpy.env.overwriteOutput = 1
infc = r'polygons'
outfc = r'divided_polygons'

totalpolycount = arcpy.GetCount_management(infc)
polygoncount = 0
partcount = 1

with arcpy.da.SearchCursor(infc,['SHAPE@AREA','SHAPE@LENGTH','SHAPE@XY','SHAPE@']) as scursor1:
    for srow1 in scursor1:
        polygoncount+=1        
        print "Processing polygon {} of {}".format(polygoncount,totalpolycount)
        #Find out if polygon is round (compactness >=0.045) or elongated
        compactness = srow1[0]/(srow1[1]**2)
        if compactness >= 0.045:
            print "It's round, dividing like a cake"
            dist=srow1[1] #/2
            radian=0.0
            trailingcoord = [srow1[2][0],srow1[2][1]+dist]            
            for j in range(1,6):
                unloopcounter=0
                areadiff=99999.0
                #Keep increasing angle of cakepiece until difference in area between piece and 1/5 of polygon is below 1/2000 of polygon area 
                while abs(areadiff)>(srow1[0]/2000) and unloopcounter<5000:
                    radian+=0.002
                    feature_info = [list(srow1[2]), trailingcoord]
                    if j<5:
                        feature_info.append([srow1[2][0]+dist*math.sin(radian),srow1[2][1]+dist*math.cos(radian)])
                    else:
                        feature_info.append([srow1[2][0],srow1[2][1]+dist])
                        unloopcounter=9999
                    polygontemp=arcpy.Polygon(arcpy.Array([arcpy.Point(*coords) for coords in feature_info]))
                    arcpy.Clip_analysis(srow1[3],polygontemp, r'in_memory\polygonpart')
                    with arcpy.da.SearchCursor(r'in_memory\polygonpart',['SHAPE@AREA']) as scursor2:
                        for srow2 in scursor2:
                            areadiff=srow1[0]/5-srow2[0]
                    unloopcounter+=1
                print 'Part {} took {} attempts'.format(j, unloopcounter)
                trailingcoord=[srow1[2][0]+dist*math.sin(radian),srow1[2][1]+dist*math.cos(radian)]
                arcpy.CopyFeatures_management(r'in_memory\polygonpart','Polygon_'+str(polygoncount)+'_Part_'+str(j))
        else:
            print "It's elongated, dividing like a bread loaf"
            arcpy.MinimumBoundingGeometry_management(srow1[3],r'in_memory\bounding','RECTANGLE_BY_WIDTH')
            arcpy.SplitLine_management(r'in_memory\bounding', r'in_memory\splitline')
            linelist=[]
            with arcpy.da.SearchCursor(r'in_memory\splitline',['SHAPE@LENGTH','SHAPE@']) as scursor3:
                for srow3 in scursor3:
                    linelist.append(srow3)
                linelist=sorted(linelist,key=lambda x: x[0])
                deltax = linelist[3][1].firstPoint.X-linelist[3][1].lastPoint.X
                deltay = linelist[3][1].firstPoint.Y-linelist[3][1].lastPoint.Y
            with arcpy.da.UpdateCursor(r'in_memory\bounding',['SHAPE@XY']) as ucursor1:
                for urow1 in ucursor1: 
                    shiftx=0.0
                    shifty=0.0
                    for i in range(1,6):
                        areadiff = 99999
                        unloopcounter = 0         
                        arcpy.Clip_analysis(srow1[3], r'in_memory\bounding', r'in_memory\boundingclip')
                        #keep moving bounding polygon until the differnce between piece left after erase and 1/5 of polygon is below 1/2000 of polygon area 
                        while abs(areadiff) > (srow1[0]/2000) and unloopcounter < 5000:
                            if i==5:
                                shiftx+=deltax
                                shifty+=deltay
                                unloopcounter=9999
                            ucursor1.updateRow([[urow1[0][0]+shiftx,urow1[0][1]+shifty]])
                            shiftx+=deltax/5000
                            shifty+=deltay/5000
                            arcpy.Erase_analysis(r'in_memory\boundingclip', r'in_memory\bounding', r'in_memory\polygonpart')
                            with arcpy.da.SearchCursor(r'in_memory\polygonpart',['SHAPE@AREA']) as scursor4:
                                    for srow4 in scursor4:
                                        areadiff = (srow1[0]/5)-srow4[0]
                            unloopcounter+=1                            
                        print "Part {} took {} attempts".format(i,unloopcounter)
                        arcpy.CopyFeatures_management(r'in_memory\polygonpart', out_feature_class='Polygon_'+str(polygoncount)+'_Part_'+str(i))
print "Merging all the parts together, deleting temp data and repairing small gaps"
partlist = arcpy.ListFeatureClasses(wild_card='Polygon_*')
arcpy.Merge_management(partlist, outfc)
for part in partlist:
    arcpy.Delete_management(part)
arcpy.Integrate_management(outfc, 0.001) #The dividing sometime creates small gaps between parts. Gaps smaller than 0.001 m are repaired (change the 0.001 if needed).
print "Finished"

Answer 4

Not sure how it works in QGIS, but in ArcGIS, there is a tool called Slice that can serve your cause. However you need to transform the file to raster to do that.

Answer 5

This could be achieved in 2 simple steps

Step 1 :

Create grid polygons:

(Use Vector>>Research Tools>>Vector grid) create a grid of polygons the same extent as your shapefile, with the right distance between divisions ('parameter')

Tip: Check with reference lines first then create polygons.

Step 2:

Intersect the two layers:

(Vector>>Geoprocessing Tools>>Intersect) first layer as the original shapefile and the second as your vector grid. The output will be your shapefile chopped up by the boundaries of the vector grid.

Thanks.