ArcGIS Desktop – Fill Polygon with Points for Maximum Density (Bin Packing Problem)

Question

Is there a tool for ArcMap 10 that fills a given polygon with point geometries?

The points should have a pre-defined distance between themselves, but the points position can vary. Fishnet is not an option because it wont create points based on the polygon's shape.

I first thought of using the Network Analyst extension but that seems to need a network layer, which doesn't exist in my case.

Here is an example of the fill-pattern I'm thinking of. Now let's say the minimum distance between each point is 100 meters. Some adjacent points have a greater distance because of the shape of the polygon

I hope it's a bit clearer now what I'm searching for

Answer 1

EXPERIMENT:

I placed points at 200m apart using script from this post and extent much bigger than polygon of interest:

I've made it topmost layer in the current mxd table of content.

I placed polygon layer below and finally created empty point feature class and made it 3rd from the top. These 3 layers are inputs to the script below.

RESULT:

Shows one of many possible solutions, where point count increased from 19 to 24:

As I mentioned in my comments, there are 3 parameters to optimise. I don’t have scipy installed this is why I applied following tactic:

1. Define near point (pClose) on the polygon outline for every point outside polygon and within 200m. Calculate coordinate shifts (dX,dY)
2. Shift all points by dX, dY
3. Rotate all new points around pClose, find angle which result in maximum point count inside polygon
4. Apply best point and angle to original dataset

One of the application is optimisation of pivot irrigation system. In this case coverage increased by 16%.

SCRIPT:

import arcpy, traceback, os, sys,math
from math import radians,sin,cos

mxd = arcpy.mapping.MapDocument("CURRENT")
layers=arcpy.mapping.ListLayers(mxd)
[triPoints,pGonLayer,outFC]=layers[:3]
L=200
gr=(math.sqrt(5)-1)/2


try:
    def showPyMessage():
        arcpy.AddMessage(str(time.ctime()) + " - " + message)

    #golden section to find minimum
    def gss(a,b,tol):
        c=b-gr*(b-a)
        d=a+gr*(b-a)
        while abs(c-d)>tol:       
            fc=f(c);fd=f(d)
            if fc<fd:
                b=d
                d=c
                c=b-gr*(b-a)
            else:
                a=c
                c=d
                d=a+gr*(b-a)
        return (b+a)/2
    # rotate points
    def rPoints(angle):
        a=radians(angle)
        movedPoints=[]
        for p in allPoints:
            x,y=p.X-pClose.X,p.Y-pClose.Y
            xN=cos(a)*x+sin(a)*y
            yN=-sin(a)*x+cos(a)*y
            pN=arcpy.Point(xN+pClose.X,yN+pClose.Y)
            if pgon.distanceTo(pN)==0:
                movedPoints.append(pN)
        return movedPoints
    #function to minimise
    def f(a):
        inside=len(rPoints(a))
        return len(allPoints)-inside

#   get polygon
    with arcpy.da.SearchCursor(pGonLayer, "SHAPE@") as rows:
        for row in rows:pgon=row[0]
#   get points inside and nearby
    pointsInside=[];nearPoints=[]
    with arcpy.da.SearchCursor(triPoints, "SHAPE@") as rows:
        for row in rows:
            p=row[0].firstPoint
            D=pgon.distanceTo(p)
            if D==0:pointsInside.append(p)
            elif D<=L:nearPoints.append(p)
        allPoints=pointsInside+nearPoints
#   iterate through near points
    pBoundary=pgon.boundary()
    nMax=len(pointsInside)
    arcpy.AddMessage("Original count of points %i" %nMax)
    for p in nearPoints:
        chainage=pBoundary.measureOnLine(p)
        pClose=pBoundary.positionAlongLine(chainage).firstPoint
        angle=gss(-60.0,60.0,0.01)
        nCur=len(rPoints(angle))
        if nCur>nMax:
            nMax=nCur;bestAngle=angle;bestPoint=pClose
    arcpy.AddMessage("\nCount of %s achieved at %i degrees angle\n" %(nMax,bestAngle))
# transfer results
    pClose=bestPoint
    movedPoints=rPoints(bestAngle)
    arcpy.AddMessage(len(movedPoints))
    curT=arcpy.da.InsertCursor(outFC,"SHAPE@")
    for p in movedPoints:
        curT.insertRow((p,))
    del curT, mxd
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()