# Finding geographic median in ArcGIS Desktop?

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.

• 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
• 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

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 = []

# 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"