# How to create quadralinear polygons around a set of points

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.

• I can see neither points or polygons. Commented Dec 21, 2018 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. Commented Dec 21, 2018 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. Commented Dec 21, 2018 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. Commented Dec 22, 2018 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. Commented Dec 22, 2018 at 4:15

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
vertDist,horMesh=getDeltas(0)
horDist,vertMesh=getDeltas(1)
buffs=[vertDist,vertDist,horDist,horDist]

arcpy.CopyFeatures_management(horMesh+vertMesh,outFile)
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! Commented Dec 24, 2018 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. Commented Dec 24, 2018 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 Commented Dec 27, 2018 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. Commented Dec 28, 2018 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! Commented Dec 28, 2018 at 5:55