shapefile - Fitting known size polygons into irregular polygons using ArcGIS Desktop?

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

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