Fitting known size polygons into irregular polygons using ArcGIS Desktop?

Question

I have a task that sounds simple, but I just can't get through the last step. I need to find suitable areas where wind farms that have a size of 500x150 meters, could be placed. The restrictions are slope - best areas are form 0-2 degrees and 2-3 degrees are worth considered - and a 3 km buffer from the selected location. Data used here was SRTM as this is preliminary study. My workflow was to buffer the locations, clip the DEM, create slope, reclass into 3 categories and convert to shapefile.

My categories of reclass are as follows:

1. First category - 1 = 0-2 degrees - green;
2. Second category - 2 = 2-3 degrees - orange;
3. Third category - 3 = all other data - red.

The issue that I am facing is that I would like know many 500x150 meters polygons can I fit in the green area, without overlapping the orange and red areas and without overlapping the polygon itself. The end result should be a shp with the possible positions for the aforementioned 500x150 meters polygon.

I should also specify that the 500x150 meters polygon can have any position/orientation and that I have access to ArcGIS Desktop 10.6/ArcGIS Pro 2.2/Global Mapper 19 and QGIS and my coding skill relate to very simple tasks that I replicate from other users.

I had a look at Checking if polygon fits inside another polygon using ArcGIS Desktop? where a similar things is described, but with circles and using ETGeowizard, also an older version of Arcmap.

Also, Packing Polygons within polygon using ArcGIS Desktop? is almost what I need, but the answers are just to above my level.

EDIT: the purple area is an existing possible location drawn by hand. I am trying to figure out the maximum number of 500x150m polygons that can be mosaic-ed into the green area. They don't necessarily have to touch each other.

Answer 1

It is very complicated task known as bin packing problem.

The script below produces one of countless sub-optimal solutions. Algorithm:

• places fish net over rotated POLYGON to find out rotation angle in range (0,175,5) that result in maximum count of complete rectangles
• breaks if no such rectangles found, otherwise
• un-rotate every good rectangle and append it to a big list
• erase original polygon by big list. POLYGON = erase result, repeat

It works with one single-part polygon, holes Ok. Results can be slightly improved, by optimizing fish net origin location, similar to this. It would result in one more cell (near 85), but I don't think it is enough value for efforts.

import arcpy
from arcpy import env
from math import radians,sin,cos
env.overwriteOutput = True
infc=arcpy.GetParameterAsText(0)
outFC=arcpy.GetParameterAsText(1)
d=arcpy.Describe(infc);SR=d.spatialReference
W=50;L=15;A=0.99*W*L
fnet="in_memory/fnet"
erased="in_memory/fnet"

# rotate polygon
def ShapeMake(pGon,angle):
    a=radians(angle)
    ARR=arcpy.Array()
    cX=cPoint.X;cY=cPoint.Y
    for part in pGon.boundary():
        ar=arcpy.Array()
        for p in part:
            x,y=p.X-cX,p.Y-cY
            xN=cos(a)*x+sin(a)*y
            yN=-sin(a)*x+cos(a)*y
            pN=arcpy.Point(xN+cX,yN+cY)
            ar.add(pN)
        ARR.add(ar)
    pgonRotated=arcpy.Polygon(ARR,SR)
    return pgonRotated
# create fishnet and count complete polygons
def fnetMake():
    FNET=[]
    ext=rotated.extent
    oc='%s %s' %(ext.XMin,ext.YMin)
    ya='%s %s' %(ext.XMin,ext.YMax)
    cc='%s %s' %(ext.XMax,ext.YMax)
    arcpy.CreateFishnet_management(fnet, oc,ya, W, L,"","",
                                   "","NO_LABELS", rotated,"POLYGON")
    rects=arcpy.Clip_analysis(fnet, rotated, g)
    for chop in rects:
        if chop.area<A:continue
        FNET.append(chop)
    return FNET

g=arcpy.Geometry()
PGON=arcpy.CopyFeatures_management(infc,g)[0]
theList=[PGON];bigList=[]

nBefore=0
while True:
    for toCut in theList:
        ## FIND rotation to maximise complete rectangles
        nMax=0
        cPoint=toCut.centroid
        for i in range(36):
            angle=5*i
            rotated=ShapeMake(toCut,angle)
            squares=fnetMake()
            N=len(squares)
            if N<=nMax:continue
            nMax=N
            keepers=squares[:]
            bestAngle=angle
        if nMax==0:continue
        arcpy.AddMessage("%s cell(s) found so far" %nMax)
        for item in keepers:
            rotated=ShapeMake(item,-bestAngle)
            bigList.append(rotated)
    if nBefore==len(bigList):break
    nBefore=len(bigList)
    arcpy.Erase_analysis(PGON, bigList, erased)
    theList=arcpy.MultipartToSinglepart_management(erased, g)
arcpy.CopyFeatures_management(bigList,outFC)