5

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.

enter image description here

  • 5
    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 6 '16 at 7:27
  • 1
    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 6 '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 6 '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 6 '16 at 15:58
  • 1
    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 6 '16 at 17:37
2

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 6 '16 at 21:31
  • Note I updated script, a bit gone missing during copy-paste. Use latest version – FelixIP Mar 6 '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 7 '16 at 15:45
  • I swapped out the pic above with a pic of the 'links' layer symbolized by individual part # – Thomas Lee Mar 7 '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 7 '16 at 21:33

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.