How to create quadralinear polygons around a set of points

Question

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:

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.

Answer 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:

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.:

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

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