# Graph/Network building and analysis of linked polygons in ArcMap

This is more of a follow-up question to Dividing polygons into *n* number of groups of equal counts with ArcPy? and Grouping village points based on distance and population size using ArcMap, particularly the answers by FelixIP.

I am trying to group a set of polygons spatially with a max population constraint. (The Grouping Analysis tool in ArcGIS is similar to desired results, but can only do temporal or spatial constraints, no population -- so some network analysis with ArcPy is needed.)

I am mostly following the steps outlined in FelixIP's answer in the first question, but am stuck at the part where he suggests "[fi] and [ti] are sequential number of connected nodes. To populate this table search this forum on how to assign from and to nodes to link." Unfortunately, there seems to be no such answer (or I'm phrasing poorly), hence my question:

How do I create a network/graph that can be used to calculate the number of sequentially connected nodes? And in fact, do I correctly understand that these [fi] and [ti] fields represent the total number of nodes that each node can be reached from (fi) and the number of neighbor nodes each node can reach (ti)?

So far I have everything else needed until that point: a Point node layer, and a Polyline link layer (with from and to node). Visual aid below for context:

(source: maks.is)

• Have a look at gis.stackexchange.com/questions/125090/… on how to assign node ids to line. It is very useful technique when working with networks. However in my solution with grouping creating lines and nodes was a massive overkill for visual presenration. In general think of your polygons as nodes, joined polygons - links. Join has all fields y9u'll need. Orig_fid = fi (from node index), etc I'll post extended answer on Monday Jul 3, 2015 at 6:31
• I'm missing what [fi] and [ti] are then. Your table had both [fi], [ti], but also [fromn] and [ton]. If all I need is which polygons each link connects, then why doesn't the [fromn] and [ton] suffice? Your example also had entries where [fi] and [ti] are the same, but you wouldn't have links connecting the same node to itself, right? Looking forward to your extended answer. Thanks! Jul 3, 2015 at 11:03

As I said, points and lines are overkill. Same with "fromn" and "ton" - from and to nodes names. This is a simplified answer, let's make it exercise.

Create shapefile like this:

Call this layer "nodes" in the table of content.

Add spatial join to itself (remove all of the fields!):

Links table will look like that

Add output to map, call it "links". Remove all the features, where "TARGET_FID"="JOIN_FID" (former fi and ti fields). Add field "times" to links table, populate it with 1. This is cost of travel from neighbour to neighbour.

Create field "P2013" (population field) in nodes layers, populate it by 1 or 2 as shown on the first picture.

Create field "rcvnode", long integer, - receiving node to store host's "FID" value.

Attach this script to custom tool:

``````import arcpy, traceback, os, sys
import itertools as itt
sys.path.append(r'C:\Users\felix_pertziger\AppData\Roaming\Python\Python27\site-packages')
import networkx as nx
RATIO = int(arcpy.GetParameterAsText(0))
tolerance=float(arcpy.GetParameterAsText(1))
try:
def showPyMessage():
arcpy.AddMessage(str(time.ctime()) + " - " + message)
# FIND LAYERS
mxd = arcpy.mapping.MapDocument("CURRENT")
theNodesLayer = arcpy.mapping.ListLayers(mxd,"nodes")[0]
G=nx.Graph()
for row in cursor:
(f,t,c)=row
pops=[row[0] for row in arcpy.da.TableToNumPyArray(theNodesLayer,("P2013"))]
length0=nx.all_pairs_dijkstra_path_length(G)

nNodes=len(pops)
aBmNodes=[]
aBig=xrange(nNodes)
host=[-1]*nNodes

while True:
RATIO+=-1
if RATIO==0:
break
aBig = set(aBig).difference(set(aBmNodes))
p=itt.combinations(aBig, 2)
nBig=len(aBig)
nTotal=nBig*(nBig-1)/2
arcpy.SetProgressor("step", "", 0, nTotal,1)
pMin=1000000
for a in p:
arcpy.SetProgressorPosition()
S0,S1=pops[a[0]],pops[a[1]]
others=set(aBig).difference(set(a))
for i in others:
p=pops[i]
L0=length0[a[0]][i]
L1=length0[a[1]][i]
if L0<L1:S0+=p
else: S1+=p
if S0!=0 and S1!=0:
sMin=min(S0,S1)
sMax=max(S0,S1)
df=abs(float(sMax)/sMin-RATIO)
if df<pMin:
pMin=df
aBest=a[:]
arcpy.AddMessage('%s %i : %i = %.2f' %(aBest,sMax,sMin, float(sMax)/sMin))
if df<=tolerance:break
lSmall,lBig,S0,S1=[],[],0,0
for i in aBig:
p=pops[i]
L0=length0[aBest[0]][i]
L1=length0[aBest[1]][i]
if L0<L1:
lSmall.append(i)
S0+=p
else:
lBig.append(i)
S1+=p

aBmNodes,m,n=lSmall[:],0,1
if S0>=S1:
aBmNodes,m,n,lBig = lBig[:],1,0,lSmall[:]
for i in aBmNodes:
host[i]=aBest[m]
for i in lBig:
host[i]=aBest[n]
# save results in nodes' table
with arcpy.da.UpdateCursor(theNodesLayer, "rcvnode") as cursor:
i=0
for row in cursor:
row[0]=host[i]
cursor.updateRow(row)
i+=1
del row, cursor
except:
message = "\n*** PYTHON ERRORS *** "; showPyMessage()
message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()

``````

Tool has 2 parameters:

Sub-optimal, one of many, local or whatever solution:

I hope this will help better understand meaning of fields in above approach

• To help script to find solution quickly, reorder your nodes so southern most has "FID"=0 Jul 5, 2015 at 23:44
• Fantastic! With tweaking to suit my data and some nuances, this worked. Albeit, with 350 parcels and 39 desired clusters, low tolerance values do not bode well. With sorting and tweaking it should do a good job. Just to know: I see that the distances algorithm is Djikstra's, but is the splitting part also a particular one, i.e. can I find info about it online to tweak for my data? Jul 6, 2015 at 11:57
• There are some rules to follow on this site: 1. If problem solved MARK it as solved. I feel sick when this is not a case 2. If you want to attract attention add @ before his/her nickname, e.g @Emil Brundage. The split technique I've cooked myself, chances are someone else invented it earlier. The best name for it is bipolar dissect. Jul 6, 2015 at 20:04
• Cheers. Marked as solved. Jul 7, 2015 at 8:09