ArcMap fill polygon with points (highest possible number) (Bin Packing Problem: Best Fit)

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

Answer 2

Since you asked specifically if there is a way to do this with out-of-the-box ArcGIS tools, here is a solution to do it with existing tools. This could also be automated with scripting.

This is the workflow I used.

1. Use Multiple Ring Buffer tool, creating negative buffers, say -100, -200, -300, etc. until you cover breadth of polygon.
2. Use Feature To Line tool to convert polygons to lines, with original polygon and negative buffers as input. Include both the original polygon and the buffers as input.
3. Project the resulting line feature class, if necessary. You want the horizontal units to be meters or feet.
4. Create an Id field (Long) and populate it (just copy ObjectId).
5. Use Create Routes tool (use distance as the measure).
6. Create an event table, creating points at the desired distance along the routes.
7. Use the Make Route Event Layer tool to create points based on your event table. Make sure you check the box to create an error column.
8. Select the points you want to keep , LOC_ERROR = NO ERROR. You may need to remove points too close together at the end/start of the lines.
9. Use the Feature To Point tool, with the inside-most buffer as the input, to create a point in the space left in the inner-most buffer. The result should be a single point in the center of the original polygon.
10. Merge the centroid and the event points into one feature class, using the Merge tool.

This is what the result looked like. I also showed the negative buffers, just to show better how it all fits together.

Answer 3

How about the Create Random Points tool (Data Management Tools > Feature Class)?

http://resources.arcgis.com/en/help/main/10.2/#/Create_Random_Points/00170000002r000000/