5

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.

Example

  • 2
    Two points to clarify: 1. The purple area is an existing wind farm so do not overlap it or the unsuitable areas in orange and red. Is that interpretation correct? 2. Are you trying to figure out the maximum number of 500x150m polygons that can be mosaic-ed into the green area or are you trying to figure out all the different possible locations for a single 500x150m within the green area? – Hayden Elza Sep 22 '18 at 0:20
  • @FelixIP If I use your suggestion I will get only one polygon per area, at least this is how I understand it, because my green area is just one large polygon and you are using centroids to determine where the polygon should be placed. – alecsx Sep 22 '18 at 6:46
  • Yes, it's a case unfortunately. This is no simple task at all. I'd test few remotest points inside polygons, if you need to find certain number of them. Hopefully. – FelixIP Sep 22 '18 at 7:00
  • As a further constraint, the centre point of one of your rectangles cant have a red/yellow pixel within 75m (half the width) of it, so if you make a 75m raster buffer on the red/yellow you'll get a smaller set of possible centre location pixels. That might be useful for an algorithm which starts by selecting possible centre locations. – Spacedman Sep 22 '18 at 7:43
7

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.

enter image description here

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.