How to group adjacent polygons with similar orientations?

Question

I have a layer of about 40k polygons, each representing an individual agricultural field. I want to group adjacent fields (within 20meters) that are within 15degrees of orientation of each other. I need to be able to measure how many of the 40k fields can be placed into groups with other fields that are contiguous and parallel. I have already calculated the orientation of each field and converted them to a range of 0-180. I am fairly familiar with ArcMap, however, I am just recently learning to apply my limited coding skills.

How would you want to handle a situation in which three fields are lined up; the western field has an orientation of 0 degrees, the center field has an orientation of 14 degrees and the eastern field has an orientation of 28 degrees? Would they all three be joined, even though two of them differ by 28 degrees? If not, how would you determine which of the outside fields the center field pairs with? — Tom, Mar 06 '16 at 07:27
In the image above on the eastern side what was the rule that grouped the almost horizontal 5 green fields as the top and bottom are not adjacent, also there are 2 fields that are twice the length. If area is not a factor then why immediately above these are the brown fields not part of the blue as they are all in the same direction. There are rules you have not explained. — Hornbydd, Mar 06 '16 at 12:23
I think a situation with multiple fields of increasing difference in orientation would be acceptable as long as the immediately adjacent fields are within the range. So these three fields you describe could go into a single field group. The picture is only symbolized by the Orientation column within the layer. It is just there to give a visual of what the layer looks like. Fields that look like they should be grouped and are not probably fall just outside the range for that color setting. — Thomas Lee, Mar 06 '16 at 15:52
Also, end-to-end fields should not be included, however, in some instances less than 50% of one field may be lined up parallel to another field of a much larger size. I do not have any rules for excluding fields that are not lined up exactly adjacent or are not of equal size. — Thomas Lee, Mar 06 '16 at 15:58
The fields are Pre-Columbian archaeological earthworks with an unknown chronology. The greatest association that can be made beyond an individual field is between two immediately adjacent fields. Adjacent meaning close enough to prevent another field from having been built between the fields and therefor preventing a farmer from expanding upon a particular field area at a later date. The groups of parallel adjacent fields can potentially represent social institutions controlling construction of the fields. So my rules are bendable and can be tweaked after I figure out how to implement them — Thomas Lee, Mar 06 '16 at 16:06
I'm curious to know how you defined orientation, since they're not all rectangular; i.e., some of them bend. The logic you're already using could be helpful in figuring out how to identify those that are side-by-side vs end-to-end. Also, are your polys densified (with vertices along straight edges, not marking a corner/bend), or are the vertices only at the corners and the few legit bends? — Tom, Mar 06 '16 at 17:35
Sorry to bombard you with so many questions, but with such an amorphous problem/approach, we've got to find some starting point. — Tom, Mar 06 '16 at 17:37

FelixIP · Accepted Answer · 2016-03-07T19:24:51.803

I suggest this workflow:

Spatially join your polygons to itself, one to many, 20m distance. In joined table create field "same", type short and compare directions using this field calculator: expression:

deg=math.pi/180
def similar(a,b,da):
 L=[math.sin((a-da)*deg),math.sin((a+da) *deg)]
 L.sort()
 return (math.sin(b*deg)>=L[0])*(math.sin(b*deg)<=L[1])

'============================

similar( !A!, !B!,15)

Delete records with !Same!=0 Call spatial join "links" and run this script from mxd:

import arcpy, traceback, os, sys, time
import networkx as nx

try:
    def showPyMessage():
            arcpy.AddMessage(str(time.ctime()) + " - " + message)
# FIND ENVIRONMENT TABLE
    mxd = arcpy.mapping.MapDocument("current")
#  GET LINKS LAYER AND FROM and TO indexes
    theLinksLayer = arcpy.mapping.ListLayers(mxd,"links")[0]
    fldFROM="TARGET_FID"
    fldTO="JOIN_FID"
#  CREATE GROUPS
    G=nx.Graph()
    arcpy.SelectLayerByAttribute_management(theLinksLayer, "CLEAR_SELECTION")
    cursor= arcpy.da.TableToNumPyArray(theLinksLayer, (fldFROM,fldTO))
    for f,t in cursor:
        G.add_edge(int(f),int(t))
    del cursor

    dictFeatures = {}; m=0
    while True:
        s=G.nodes()[0]
        A=list(nx.ancestors(G,s))
        A.append(s)
        for i in A:dictFeatures[i]=m
        G.remove_nodes_from(A)
        n=G.number_of_nodes()
        if n==0: break
        m+=1
# TRANSFER DATA TO TABLE
    try:
            arcpy.AddField_management(theLinksLayer, "PART", "SHORT")
    except:pass
    arcpy.AddMessage('\nTransfering results to table\n')
    with arcpy.da.UpdateCursor(theLinksLayer,(fldFROM,"PART")) as cursor:
        for fid,prt in cursor:
            prt=dictFeatures[fid]
            cursor.updateRow((fid,prt))

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()

Script will create field "PART" to store group numbers, you'll need to transfer this info back to "nodes", i.e. name your polygons nicely beforehand.

Script using module networkx, let me know if installation of networkx is going to be an issue, I'll dig out script which does it without networkx

UPDATE I overcomplicated things with orientation trying to handle situation 0/180 degrees which are the same to me. Unfortunately it makes 45/135 also the same, which is not the case. Use this to find similar orientation links:

def similar(a,b,da):
 L=[a,b]
 return (max(L)-min(L)<=da)

and I leave it to you to handle cases 0/180 or similar.

To transfer !part! field back to polygons use common fields in polygons and links, e.g. names or even FID and JOIN_FID (ADD JOIN, field calculator). Also note script tested on shapefiles.

Thanks. I am not familiar with networkx but I am going to attempt to install it now. — Thomas Lee, Mar 06 '16 at 21:31
Note I updated script, a bit gone missing during copy-paste. Use latest version — FelixIP, Mar 06 '16 at 22:09
I have completed everything up to the last step where i 'need to transfer this info back to nodes. I have a field 'PART' with about 11k group numbers, but I am still trying to figure out what you meant by the last part. I can see that in the code each polygon was passed through a loop to remove the nodes. Looking at the 'links' layer after the script has run, many of the groups have been identified, however, there are some exceptionally large groups which failed to separate out groups with greatly different orientations. — Thomas Lee, Mar 07 '16 at 15:45
I swapped out the pic above with a pic of the 'links' layer symbolized by individual part # — Thomas Lee, Mar 07 '16 at 15:52
I don't see anywhere that this answer addresses end-to-end vs side-by-side polygons. That's the part of this problem that is tripping me up. — Tom, Mar 07 '16 at 21:33
It seems he doesn't care about side-by-side etc. What is critical here is definition of orientation — FelixIP, Mar 07 '16 at 21:37
@FelixIP, from his comment: "Also, end-to-end fields should not be included" — Tom, Mar 07 '16 at 21:42
I did say that end-to-end fields should be excluded, however, getting the groups defined by proximity and orientation is more important in the end. — Thomas Lee, Mar 08 '16 at 00:02
Works like a charm! to Thanks FelixIP and everyone else who took interest! This is going to let me get to work and I can figure out how to single out end-to-end fields later on! — Thomas Lee, Mar 08 '16 at 06:44

How to group adjacent polygons with similar orientations?

1 Answers1

Linked