I have a polygon layer of watersheds (N = 1950) for South Africa. Within each watershed is a value (number of people; POPULATION field). I need to find a way to automatically group contiguous watersheds so that when a group of watersheds reaches a specific threshold (say 100 people), that group is assigned a unique code in a new attribute field called GROUP.

The code that I'm looking for would need to run through the entire polygon layer, essentially making new 'larger' groups. I'd image that in a polygon layer of 1950 watersheds, there would be ~300 unique groups.

Exactly the same question has been posted previously (Clustering polygons by sum of values using ArcPy) but unfortunately they never posted their code.