3

Is it possible to split the total population in 501 polygons, each polygon with a different number of inhabitants into two North/South sub-areas with roughly the same number of inhabitants in each sub-area?

The split line, trending east-west should be along administrative boundaries of the polygons (it is called geographic median).i.e. no polygon divided/shared between sub-areas.

The picture shows the region with 501 villages (polygons) with different numbers of inhabitants each.

enter image description here

  • 1
    Could you elaborate on your first point? It's not clear what you want to do (expected output) or why? – GISKid Apr 4 '18 at 19:17
  • I would like to split the region ( that you see in the picture ) in two areas, with the same number of inhabitants in both areas. But I need to keep it within the administrative boundaries of the villages. – Aneta Malafová Apr 4 '18 at 19:31
  • 1
  • Thank you! It looks good but It is too complicated for me so I don't understand what should I do exactly. I need step by step. LOL – Aneta Malafová Apr 4 '18 at 20:35
  • See if this helps gis.stackexchange.com/questions/153094/… – FelixIP Apr 5 '18 at 0:20
4

In attempting to develop a solution to Geographic median in ArcGIS, I came up with the script below. May not be sophisticated or pythonic but seems to get what I want.

The steps:

  1. Take the shapefile(1) in question and add a centroid_(Y or X) field depending on desired direction of split
  2. Calculate total_pop and half of the total_pop (this is the marker to stop the for loop). Assumes pop field exits.
  3. Access the shapefile(1) again and create a new shapefile(2) with sorted centroid values (ascending or descending)
  4. use the new sorted shapefile(2) to find the centroid value (median value) of the feature that tips the progressively summed pop field of features to over half
  5. select all features with centroids less than this value(median value) and create another new shapefile(3).

The total pop of features in this last(half_features.shp) shapefile(3) adds up to about half the total in the original shapefile(1).

My understanding of the question is that the OP was looking for median line as shown here, except with a line running along the boundaries.

    import arcpy, os

    path = r"C:\folder\folder2"

    arcpy.env.workspace = path
    arcpy.env.overwriteOutput = True

    inputFc = "test.shp" # replace with the right file
    outFeature =  os.path.splitext(inputFc)[0] + "_lower_half.shp"

    # use CENTROID_X for east/west split
    # use CENTROID_Y for north/south split
    # POP_FIELD is the field(population) whose median line is of interest

    fields = ["CENTROID_X", "POP_FIELD"] 
    popList = []
    popList2 = []

    # Add centroid field 
    arcpy.AddGeometryAttributes_management(inputFc , "CENTROID")

    # sum pop field to get pop_total/2
    with arcpy.da.SearchCursor(inputFc, fields) as SR:
        for SearchRow in SR:
            popList.append(SearchRow[1])

    half_Pop = (sum(popList))//2

    #sort centroids (median is calculated on sorted list)
    arcpy.Sort_management(inputFc, "in_memory" + "\\" + "sortedFc", [[fields[0], "ASCENDING"]])

    # calculate median centroid value where pop sum just goes over half
    with arcpy.da.SearchCursor("in_memory" + "\\" + "sortedFc", fields) as SR:
        for SearchRow in SR:
            popList2.append(SearchRow[1])
            if sum(popList2) > half_Pop:
                continue
            medianCentroidVal = SearchRow[0]


    where_query = """ "{0}" < {1}""".format(fields[0], medianCentroidVal)

    tempLyr = arcpy.MakeFeatureLayer_management("in_memory" + "\\" + "sortedFc" , "lyr")
    arcpy.SelectLayerByAttribute_management(tempLyr, "NEW_SELECTION", where_query)
    arcpy.CopyFeatures_management(tempLyr, outFeature)

    # clean up 
    arcpy.Delete_management(tempLyr)

    mxd = arcpy.mapping.MapDocument("CURRENT")
    for df in arcpy.mapping.ListDataFrames(mxd):
        for lyr in arcpy.mapping.ListLayers(mxd, "", df):
            if lyr.name == "sortedFc":
                arcpy.mapping.RemoveLayer(df, lyr)

    del mxd

Here as an image to go along with the narrative. EDIT: some improvements based on new "learnings"

enter image description here

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.