Checking if polygon fits inside another polygon using ArcGIS Desktop?

Question

How do I check if one polygon fits inside another polygon? Furthermore how do I automatically check if one polygon fits in other polygons of my shapefile?

I have a shapefile with thousands of differnt polygons, each one listed in the attribution list with its own ID (they are singleparts, not multipart). The polygons have a non-regualar shape. Now, I want to check in which of those polygones fits one specific other polygone (it is a circle with a specific diametre). I only have a need for those polygones inside my shapefile in which this circle fits. Background (real world use): I am searching for (geothermal) sites at which I want to install a rig (with safety distances). So if there are areas where I can't place my drill rig incl. the safety distance I have no use for this area.

An idea was to have use the geometrical midpoint, overlay the circles and polygones midpoind and see if they edges overlay. Wrong idea since the geometrical midpoit is the wrong approch already. The polygons have very weird shapes, can be very long with a huge tip or whatever. Here the geometrical midpoint would be somewhere inside or even outside the long part of the polygone, the circle wouldn't fit, but at the tip it actually would fit.

Any idea how I can process automatically? Basically I just need to sort out those polygons in which my circle/other polygone doesn't fit. Don't care if the approch is for vector or raster (vector would be nice but I can always transform), for ArcGIS Desktop.

Illustration as below: In which polygons doesn't the circle fit?

Answer 1

I am using vector and raster approach to solve it. The script below has 3 input parameters: a) polygons (must be single part), b)output point feature class to store the points most distant from polygon boundary c) tolerance in map units.

Algorithm:

POLYGON =shape geometry

1. calculate centroid of POLYGON (p1)
2. Define distance from centre to polygon outline
3. Buffer outline by above distance and subtract it from polygon
4. Calculate centre of largest remaining polygon (p2)
5. Break (p2 is our point, distance is our radius) if distance between p1 and p2 <= tolerance
6. POLYGON=polygon. Go to 1

The idea borrowed if memory serves from very old post by @whuber.

import arcpy, traceback, os, sys    
from arcpy import env

env.overwriteOutput = True    
infc = arcpy.GetParameterAsText(0)    
outfc=arcpy.GetParameterAsText(1)    
tolerance=float(arcpy.GetParameterAsText(2))

d=arcpy.Describe(infc)    
SR=d.spatialReference

try:
    def showPyMessage():

        arcpy.AddMessage(str(time.ctime()) + " - " + message)

    theFolder,theFile=os.path.split(outfc)
    arcpy.CreateFeatureclass_management(theFolder, theFile, "POINT", infc, "DISABLED", "DISABLED", SR)
    arcpy.AddField_management(outfc, "theDist", "DOUBLE")
    shapefield,fidName = d.ShapeFieldName,d.OIDFieldName
    fileds2add=[]
    fields = [f.name for f in arcpy.ListFields(infc)]
    for f in fields:
        if f not in (shapefield,fidName,"Shape_Length","Shape_Area"):
            fileds2add.append(f)
    tbl=arcpy.da.TableToNumPyArray(infc,fileds2add)
    nF=len(tbl)

    arcpy.SetProgressor("step", "", 0, nF,1)
    arcpy.AddMessage("Creating points... ")
    fileds2add=["SHAPE@","theDist"]+fileds2add
    curT = arcpy.da.InsertCursor(outfc,fileds2add)
    theM=0
    with arcpy.da.SearchCursor(infc, "SHAPE@") as rows:
        for row in rows:
            feat = row[0]; prt=feat.getPart(0)
            feat=arcpy.Polygon(prt,SR)
            outLine=feat.boundary(); pSave=feat.trueCentroid
            d=outLine.distanceTo(pSave)
            if d<=tolerance:break
            while (True):
                theLine=feat.boundary()
                p=feat.labelPoint                
                d=theLine.distanceTo(p)
                try:
                    buf=theLine.buffer(d)
                except:
                    pSave=feat.labelPoint
                    break
                intR=feat.difference(buf)
                n=intR.partCount;aMax=0
                for i in xrange (n):
                    prt=intR.getPart(i)
                    pgon=arcpy.Polygon(prt,SR)
                    aCur=pgon.area
                    if aCur>aMax:
                        aMax=aCur;feat=pgon
                pN=feat.labelPoint
                d=arcpy.PointGeometry(p).distanceTo(pN)
                if d<=tolerance:
                    pSave=pN; break
            d=outLine.distanceTo(pSave)
            theRow=[pSave,d]; theP=list(tbl[theM])
            theRow+=theP
            curT.insertRow(theRow)
            theM+=1
            arcpy.SetProgressorPosition()
    del tbl

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

Script adds field theDist to point table as well as all visible attributes of polygon input feature class. Use this distance to make relevant buffer around points, or just select the points that have theDist<=your threshold.

Note it is very slow (10 polygons like one on the picture per second). Try on selection and use reasonable tolerance.

There is faster raster solution:

1. Convert polygons to outlines
2. Calculate Euclidean distance
3. Convert polygons to raster
4. Calculate Zonal statistics (max)
5. Raster calculator Con (abs(distance-max)
6. Convert to points
7. Spatial join to polygons and remove duplicates of points. Raster value field in points is radius of max circle to fit into your polygon. Control precision by cell size, does not work for overlapping polygons

UPDATE Nov.2016 The below script is modification of original. It handles donut polygons and works slightly faster, approximately 25 polygons/second:

import arcpy, traceback, os, sys
from arcpy import env
env.overwriteOutput = True
infc = arcpy.GetParameterAsText(0)
outfc=arcpy.GetParameterAsText(1)
tolerance=float(arcpy.GetParameterAsText(2))
d=arcpy.Describe(infc)
SR=d.spatialReference
try:
    def showPyMessage():
        arcpy.AddMessage(str(time.ctime()) + " - " + message)
    theFolder,theFile=os.path.split(outfc)
    arcpy.CreateFeatureclass_management(theFolder, theFile, "POINT", infc, "DISABLED", "DISABLED", SR)
    arcpy.AddField_management(outfc, "theDist", "DOUBLE")
    shapefield,fidName = d.ShapeFieldName,d.OIDFieldName
    fileds2add=[]
    fields = [f.name for f in arcpy.ListFields(infc)]
    for f in fields:
        if f not in (shapefield,fidName,"Shape_Length","Shape_Area"):
            fileds2add.append(f)
    tbl=arcpy.da.TableToNumPyArray(infc,fileds2add)
    nF=len(tbl)

    arcpy.SetProgressor("step", "", 0, nF,1)
    arcpy.AddMessage("Creating points... ")
    fileds2add=["SHAPE@","theDist"]+fileds2add
    curT = arcpy.da.InsertCursor(outfc,fileds2add)
    with arcpy.da.SearchCursor(infc, "SHAPE@") as rows:
        for m,row in enumerate(rows):
            feat = row[0]
            outLine=feat.boundary()
            pSave=feat.labelPoint
            d=outLine.distanceTo(pSave)
            while (True):
                theLine=feat.boundary()
                p=feat.labelPoint                
                dist=theLine.distanceTo(p)
                try:
                    buf=feat.buffer(-dist)
                except:
                    pSave=feat.labelPoint
                    break
                try: pN=buf.labelPoint
                except:pN=feat.labelPoint
                d=arcpy.PointGeometry(p).distanceTo(pN)
                if d<=tolerance:
                    pSave=pN; break
                feat=buf
            d=outLine.distanceTo(pSave)
            theRow=[pSave,d]; theP=list(tbl[m])
            theRow+=theP
            curT.insertRow(theRow)
            arcpy.SetProgressorPosition()
    del tbl
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

Unlike ACCEPTED ANSWER this process considers HOLES inside the polygons. e.g. If you use accepted answer then there will be error where holes exist as below image!

---

If i got you, you want to check if your circle can be inscribed inside the polygons.You can do it in ETGEOWIZARD,an arcmap addin, (see pics below) for those polygon and separate those maximum inscribed circles have area greater than or equal to that of your supplied circle-you can use select by attribute to separate those circles(>= your sample circle area) , generate centroid of those selected circles and run select by location to select those polygon that intersects with previously selected circles-center.After all it runs fast.

N.B. I am not affiliated to any software

Answer 3

In ArcGIS:

• Union the two feature classes (FCs), lets say FC1 is the main FC and FC2 is the FC you're testing to see whether it fits.
• Calculate summary statistics with statistics field being the min of the FID of FC1 and the case field being the FID of FC2.
• Join the resulting table back to FC2.
• Select by attributes where FREQUENCY = 1 and MIN FID of FC1 > -1.

Answer 4

Here is a solution on QGIS(if i got you correctly)

You should first do a copy your layer,
then install the plugin called spatial request : docs

You can now do the usual GIS Stuff ie :
select features that are completely whtin another feature, ect... see the post from Slslam !

This plugin will select these features, if you want to save them in another layer, use saveonly selected layers !

Hope that helps