Automated merge of polygons up to a limit

I have a polygon feature class with a field, let's say it's population. I need to find an automated way to merge these polygons up to a given population number.

E.g. I have polygons with population from 0 to 211. And I need to create districts with about 400 habitants. (Between 380 and 420.) I work with ArcGIS 10.

• Hmmmmmm. How should the method choose whether to merge the 211 polygon with the neighboring 170 (as opposed to choosing the 130 and 70, or 159 and 0)? Aug 7 '14 at 12:16
• I think I will need a script for this problem. As there's no if statement in ArcGIS script is needed anyway to check if a feature's population is above or below the limit. My thougth was to loop through every feature and find the neighbouring features for each and loop through the neighbours as well starting from the highest population and merging them together if the total population of the two polygons doesn't go above the limit. If it does then the script would skip that neighbour and go to the next one. Aug 7 '14 at 13:50
• You might be able to get some inspiration from this earlier question and answer of mine: gis.stackexchange.com/questions/25825/… Aug 7 '14 at 13:59
• Thank you nmpeterson! This was a great help, I was able to write a script using pysal's maxp function. It's not finished yet but I'm on the scent now! :) Aug 25 '14 at 8:19

After several tries I had a solution for the problem. I'm still not perfectly satisfied with the result but it does what it has to. It creates a list of IDs and REGION IDs. If I Join this CSV to the shapefile I can use Dissolve by REGION ID and that's it! I've created a tool of it in ArcGIS. It may help someone so I copy my code here:

import arcpy
import random
import pysal
import scipy
import numpy as np

# PARAMETERS

# The minimum value the merged polygons have to reach.
# Parameter in script: Minimum value:
# Data type: Long
FLOOR = int(arcpy.GetParameterAsText(0))

# The path of the shapefile to be merged.
# Parameter in script: Path of input shapefile:
# Data type: Shapefile
shapefile = arcpy.GetParameterAsText(1)

# The name of the ID field in the shapefile to be merged.
# Parameter in script: ID field in input shapefile:
# Data type: Field
# Obtained from: Path of input shapefile:
IDfield = arcpy.GetParameterAsText(2)

# The name of the field containg the values to sum up to FLOOR.
# Parameter in script: Field containing values to sum:
# Data type: Field
# Obtained from: Path of input shapefile:
Vfield = arcpy.GetParameterAsText(3)

# The name of the output CSV file containing polygon IDs and region codes for each input polygon.
# Parameter in script: CSV file containing region codes:
# Data type: File
# Direction: Output
# Filter: csv;
csv = arcpy.GetParameterAsText(4)

# Seed value for random.
# Parameter in script: Seed:
# Data type: Long
seed = int(arcpy.GetParameterAsText(5))

'A klaszterizálás során használt kezdő érték.'
# Initial value for clustering.
# Parameter in script: Initial:
# Data type: Long
initial = int(arcpy.GetParameterAsText(6))

# PREPARATION

# Randomization.
random.seed(seed)
np.random.seed(seed)

# Binary matrix for rook neighbourhood.
weight = pysal.rook_from_shapefile(shapefile,idVariable=IDfield)

dbf = pysal.open(shapefile[:-4] + '.dbf')
IDarray = np.array([dbf.by_col[IDfield]])
IDarray = IDarray.transpose()
Varray = np.array([dbf.by_col[Vfield]])
Varray = Varray.transpose()

# REGIONALIZATION

solution = pysal.region.Maxp(weight,Varray,FLOOR,Varray)
regions = solution.regions

# SAVE RESULT - .CSV

csv_write = open(csv,'w')
csv_write.write('"IDfield","REGION_ID"\n')
for i in range(len(regions)):
for geo in regions[i]:
csv_write.write('"' + str(geo) + '","' + str(i+1) + '"\n')
csv_write.close()