6

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:

original map split by regions

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

some regions with splitted areas

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.

  • 3
    Could you add a screenshot or better yet, upload your shapefile to a public server so that others could test it? – Joseph Oct 7 '16 at 12:17
  • 1
    I guess it's best if you draw a picture of what you have and what you want. It is --- to me --- quite unclear from the question. – til_b Oct 7 '16 at 12:21
  • 1
    In what way do you want to divide up the city into 5 regions? – Simbamangu Oct 7 '16 at 13:01
  • I'm sorry, I added some pics to better explain it. I hope it's clear now. – Rod0n Oct 7 '16 at 13:27
  • Parcel Fabric offers a Polygon, divide by equal-area. help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/… – klewis Oct 7 '16 at 14:38
7
+50

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.


  1. 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:

    Example

    And here is running the plugin:

    Running buffer 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):

    Results of plugin


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

    Model

    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):

    Model parameters


  1. 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...

    Final results

  • 1
    thanks for your answer. The problem with the solution is that the areas are one inside another and it'll be a problem for the application I'd like todo. I'm looking for something like the image I posted in the original question. This is something similar of what I want to do: bl.ocks.org/mbostock/8814734 – Rod0n Oct 7 '16 at 15:32
  • Yes, I agree that it's not as exact as you would have wanted. I would also be very interested in seeing how to achieve the same results as you wanted in QGIS =) – Joseph Oct 10 '16 at 9:37
4

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

2

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: enter image description here

If it's not round, it's elongated and will be divided like this: enter image description here

But not like this (divisions not parallel): enter image description here

Example: enter image description here

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"
  • 1
    +1, great answer even if it's for ArcMap! – Joseph Oct 21 '16 at 12:48
  • @BERA this solution is not working at and giving error at 'arcpy.SplitLine_management..' any solution ? – apaleja Apr 9 '17 at 15:48
  • What error msg do you get? – BERA Apr 11 '17 at 5:07
1

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.

  • In the Slice Documentation, Pixels of the same value are not contiguous in the reclassified raster, so would not form a single polygon. Just an observation. – klewis Oct 7 '16 at 14:33
  • Thanks. Never noticed that. Does it matter if I change the slice method? – AndrewLebron Oct 7 '16 at 15:05
1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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