Finding geographic median in ArcGIS Desktop?

Question

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.

Answer 1

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"