6

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?

overview!

In which polygons does the top left circle fit

Some fit some don't

  • Are these comparable two polygons are single part or multipart, and they in projected coordinate systems, have you tried topology? – SIslam May 19 '15 at 13:22
  • They are single part (there used to be some multipart which I transformed into single parts). I use ETRS_1989_UTM_Zone_32N so it is a projected coord. system. What do you mean by if I tried topology? Which specific toolset are you talking about? – GeoDH May 19 '15 at 13:28
  • Can u add some picture of what you need from what..i.e image of input and output.Thanks. – SIslam May 19 '15 at 13:36
  • Here is a link to a picture which explains my request, the question is in which polygons (there are many more) this red circle (top left)fits: Picture Click Here (OT: How do I post this graphic visually here?) – GeoDH May 19 '15 at 13:51
  • Ok try topological Rule "Must cover each other". But you will get what exactly fits to circle -no approximation is allowed. If u need approximation then there be another way. – SIslam May 19 '15 at 14:00
8

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.

enter image description here

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()
  • Con (abs(distance-max)<=cellSize, distance) – FelixIP May 21 '15 at 2:25
  • Since I never worked with script code inside ArcGIS I try to do this manually. I always just used the toolbox, modelbuilder, and field calculator. – GeoDH May 21 '15 at 13:27
  • This looks exactly what I am looking for btw. Thx. Just try to implement the "faster raster solution" into the modelbuilder. – GeoDH May 21 '15 at 14:11
  • To implement raster approach you don't need model builder. You can do it on entire set, unless it result in a very heavy raster – FelixIP May 21 '15 at 22:28
  • Worked just fine but needed a while to figure out how to implement the process manually. But thank you very much. Case closed!! – GeoDH May 22 '15 at 12:04
6

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!

error demo


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.

demo

enter image description here

N.B. I am not affiliated to any software

  • @slslam thanks for your efford. I don't want to use any other software than ArcGIS or QGIS so I followed your tip with the Minimum Bounding Geometry in ArcGIS. But I don't get how this should help cuz I need it the other way around, somehow. A circle around my polygons doesn't help me. If this algorithm would create the biggest possible circle inside a polygon than this would help since I than could sort after area (in the attribute table). – GeoDH May 20 '15 at 14:27
  • It workes just fine with ET. But it is a commercial software, though. But in Demo mode everything works. If you just need it once (like me) this is fine. – GeoDH May 22 '15 at 12:02
1

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.

enter image description here

  • With this method I need to place a circle per polygon from hand. I am looking for an automatic working version of what you presented in your picture above. – GeoDH May 21 '15 at 12:46
0

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 !

plugin gui

Hope that helps

  • QGIS is fine. Looks promissing. Tried it and didn't give me any results so far but I keep on trying and see if it works. – GeoDH May 20 '15 at 15:01
  • you can also use the topology tool : see here gis.stackexchange.com/questions/34415/… – julsbreakdown May 20 '15 at 15:02
  • As I mentioned in a comment on the question this would be another way to test if existing shapes fell within (it's basically the same check that the topology tool mentioned in the original answer runs), but if you're just trying to select shapes that a theoretical/representative circle will fit within, this won't work. This wasn't entirely clear in the original question, but the updated graphics better illustrate the challenge. – Chris W May 21 '15 at 0:17
  • didn't work since it is the same thing as topology tool from ArcGIS just for QGIS. – GeoDH May 21 '15 at 12:43

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.