arcgis desktop - How to group adjacent polygons with similar orientations?

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.

Answer

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.