2

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 ArcGIS Desktop?, 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:

1
(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 – FelixIP Jul 3 '15 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! – Max K Jul 3 '15 at 11:03
8

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:

enter image description here

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

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

enter image description here

Links table will look like that

enter image description here

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]
        theLinksLayer = arcpy.mapping.ListLayers(mxd,"links")[0]
        arcpy.SelectLayerByAttribute_management(theLinksLayer, "CLEAR_SELECTION")        
        linksFromI="TARGET_FID"
        linksToI="JOIN_FID"
        G=nx.Graph()
        arcpy.AddMessage("Adding links to graph")
        with arcpy.da.SearchCursor(theLinksLayer, (linksFromI,linksToI,"Times")) as cursor:
                for row in cursor:
                    (f,t,c)=row
                    G.add_edge(f,t,weight=c)
        pops=[row[0] for row in arcpy.da.TableToNumPyArray(theNodesLayer,("P2013"))]
        arcpy.AddMessage("Computing distance matrix")
        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
                arcpy.AddMessage ('Target ratio %i' %RATIO)
                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:

enter image description here

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

enter image description here

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 – FelixIP Jul 5 '15 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? – Max K Jul 6 '15 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. – FelixIP Jul 6 '15 at 20:04
  • Cheers. Marked as solved. – Max K Jul 7 '15 at 8:09
  • Sometimes ArcGis can't see layers. Select all nodes and try again. – FelixIP Jul 10 at 22:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.