0

I have a point shapefile of locations of trees in an orchard that are almost, but not quite, on a grid. The orchard is divided into plots made up of groups of trees that either in a rectangle of rows and columns or as a series of 5 trees in a single row.

I would like to create a polygon coverage of the plots where the edges of the plot are halfway between each row of trees and halfway between the last tree of one plot and the first tree of the next, so that the entire area is filled with no gaps. At the edges, the plots should extend 1/2 the row spacing and 1/2 the tree spacing beyond the trees within the plot. Essentially a body-centered lattice, like:

enter image description here

Not all plots have the same number of trees along the rows. The point data are just the location of trees. I need to create the polygons and then attribute them with the data for each plot. Each set of rows has the same number of trees and number of rows for each plot. (note the rows are defined North to South by the way the trees were planted so you could think of them as columns).

Is there a way to create these polygons without individually measuring the locations of the corner vertices? Alternatively I would settle for rectangles but that wouldn't be as rigorous.

  • 2
    I can see neither points or polygons. – FelixIP Dec 21 '18 at 4:19
  • You could have a look at Thiessen (Voronoi) polygons, perhaps this would suit your purpose but aren't guaranteed to be square. resources.arcgis.com/en/help/main/10.2/index.html#//… (advanced license required). In a fairly regular grid they will be close to square but are guaranteed to cover the entire area between the rows. – Michael Stimson Dec 21 '18 at 4:22
  • While using the fixed-width formatting of code block would make you ASCII art more legible, firing up Paint on a screenshot of your data would make your point much clearer. Don't forget to put the exact software release in the body of the question. – Vince Dec 21 '18 at 5:03
  • I have added a jpg of crude polygons, showing what I am after. My current thought is to draw the polygons by snapping to the tree points and then shift over and up 1/2 the tree and row spacing. I don't know how to automate that though. – haresfur Dec 22 '18 at 0:49
  • Now please provide the software. You've tagged this with 10.0, which has been retired for many years. It's doubtful that any modern programming with arcpy would function on 10.0, which would make this task much more difficult. – Vince Dec 22 '18 at 4:15
1

I wrote a script that estimates average distances between rows and columns of points. It creates 2 "horizontal" and 2 "vertical" sides for points' minimum bounding rectangle and populate sides table with the half of the relevant distance:

enter image description here

Orientation of rows doesn't matter, it will handle a missing tree(s). When run from mxd, script will honor selection.

I created parallels to the sides using solution from this post.:

enter image description here

and use them as input to compute fish net. Solution.:

enter image description here

SCRIPT:

import arcpy
from arcpy import env
env.overwriteOutput = True
lyr="points"
MBG="in_memory\\mbg"
outFile='...\\sides.shp'

## DRAW MESH  
def getSetOfLines(dL,first):
    A=sides[first]; B=sides[first+2]
    L2=A.length; L=dL
    horizontal=filter(lambda x: x not in [A,B],sides)
    while True:
        p1=A.positionAlongLine (L).firstPoint
        p2=B.positionAlongLine (L2-L).firstPoint
        line=arcpy.Polyline(arcpy.Array([p1,p2]),SR)
        horizontal.append(line)
        L+=dL
        if L>L2:break
    arcpy.Near_analysis(MBG, horizontal)
    distances=[row[0] for row in arcpy.da.TableToNumPyArray(MBG,"NEAR_DIST")]
    retValue=sum(distances)/len(distances)
    return (retValue,horizontal[:2])

## FIND BEST SET BY INCREASING NUMBER OF SPLITS
def getDeltas(f):
    sMax=1e12
    N=2
    lSide=sides[f].length
    while True:
        sCur=getSetOfLines(lSide/N,f)[0]
        if sCur>sMax:break
        arcpy.AddMessage("Average distance to mesh %s" %int(sCur))
        N+=1; sMax=sCur
    N-=1
    return (lSide/N/2,getSetOfLines(lSide/N,f)[1] )

# MAIN MODULE
d=arcpy.Describe(lyr)
SR=d.spatialReference
g=arcpy.Geometry()
shp=arcpy.MinimumBoundingGeometry_management(lyr,g, "RECTANGLE_BY_WIDTH", "ALL")[0]
part=shp.getPart(0)
sides=[]
for i in range(1,len(part)):
    cut=part[i-1:i+1]
    sides.append(arcpy.Polyline(cut,SR))
## re-use MBG feature class
arcpy.CopyFeatures_management(lyr,MBG)
## FIND STEPS
arcpy.AddMessage("Finding mesh dimensions")
vertDist,horMesh=getDeltas(0)
horDist,vertMesh=getDeltas(1)
buffs=[vertDist,vertDist,horDist,horDist]

arcpy.CopyFeatures_management(horMesh+vertMesh,outFile)
arcpy.AddField_management(outFile,"BUFFER","Double")
N=0
with arcpy.da.UpdateCursor(outFile,"BUFFER") as cursor:
    for row in cursor:
        row[0]=buffs[N]
        N+=1
        cursor.updateRow(row)

arcpy.RefreshActiveView()

UPDATE: try to replace relevant line by:

cut=arcpy.Array([part.getObject(i-1),part.getObject(i)])
  • That looks like what I need. provided it can be tweaked to deal with plots that consist of 3 rows of 4 trees, like the western side of my illustration. Thanks! – haresfur Dec 24 '18 at 3:55
  • Number of points inside polygon is for you to define at step 3. In my example it could be 1 - 270 tiles, 135 - 2 tiles, 270 - 1 tile. I was hoping it is clear, when described process. – FelixIP Dec 24 '18 at 5:03
  • plowing through this with my meger python. I get to for i in range(1,len(part)): cut=part[i-1:i+1] sides.append(arcpy.Polyline(cut,SR)) and it throws the error: Runtime error Traceback (most recent call last): File "<string>", line 59, in <module> File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\arcobjects\mixins.py", line 161, in getitem return convertArcObjectToPythonObject(self._arc_object.GetObject(index)) AttributeError: Array: Error in parsing arguments for GetObject – haresfur Dec 27 '18 at 23:04
  • It means that minimum bounding geometry is not a rectangle. You need at least 3*3 points selected. Also check projection to be defined. – FelixIP Dec 28 '18 at 1:23
  • Hmm, runs through the "sides=[ ]" statement and part appears to be an array of 5 points with the first and last points the same and coordinates are close to those of my corner tree points. It won't let me slice part even if I type part[0:2] into the interpreter. Makes no sense to me. Projection is defined and I even transformed it to be the same as the data frame. This is the first time I used python since getting a new computer, if that matters. Thanks for your time! – haresfur Dec 28 '18 at 5:55

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.