Circle filling a polygon

1341
1
Jump to solution
12-07-2022 07:01 AM
teakplantation
New Contributor III

For any given polygon in a feature class I would like to create the largest possible circle that fits within this polygon (see examples below); in case of the original polygon being a circle, the new circle would be identical to the original.

Is there a function in ArcGIS Pro to achieve this?

PS: the minimum bounding geometry tool can draw a minimum circle on the outside of each polygon (https://pro.arcgis.com/en/pro-app/latest/tool-reference/data-management/minimum-bounding-geometry.ht...). I just need to achieve the circle within the polygon instead.

0 Kudos
1 Solution

Accepted Solutions
JohannesLindner
MVP Frequent Contributor

AFAIK, there is no out of the box tool for that.

This thread lists a few algorithms the get an incircle: algorithm - Largest circle inside a non-convex polygon - Stack Overflow

 

I took one of the approaches:

  • tesselate the polygon with a Voronoi diagram (Thiessen polygons in ArcGIS)
  • find the node with the greates distance to the polygon's border, that's the center of the incircle
def create_incircle(polygon):
    add = arcpy.env.addOutputsToMap
    arcpy.env.addOutputsToMap = False
    # densify polygon to get rid of curves, which only have 2 vertices
    polygon = polygon.densify("DISTANCE", 10, 1)
    # extract vertices into a point fc
    vertices = arcpy.management.CreateFeatureclass("memory", "vertices", "POINT", spatial_reference=polygon.spatialReference)
    with arcpy.da.InsertCursor(vertices, ["SHAPE@"]) as cursor:
        for part in polygon:
            for v in part:
                cursor.insertRow([v])
    # create Thiessen polygons
    thiessen = arcpy.analysis.CreateThiessenPolygons(vertices, "memory/thiessen")
    # intersect these with  themselves to get candidate points
    test_points = arcpy.analysis.Intersect([thiessen], "memory/TestPoints", output_type="POINT")
    # get the point inside the polygon with the greatest distance to the polygon's boundary
    candidates = []
    boundary = polygon.boundary()
    with arcpy.da.SearchCursor(test_points, ["SHAPE@"]) as cursor:
        for p, in cursor:
            if p.disjoint(polygon):
                continue
            candidates.append([p, p.distanceTo(boundary)])
    center = sorted(candidates, key=lambda c: c[1])[-1]
    # buffer that point ith the distance to the boundary, return
    arcpy.env.addOutputsToMap = add
    return center[0].buffer(center[1])



# create the ouput fc
incircles = arcpy.management.CreateFeatureclass("memory", "Incircles")
arcpy.management.AddField(incircles, "FID", "LONG")

# create the incircles
with arcpy.da.InsertCursor(incircles, ["SHAPE@", "FID"]) as i_cursor:
    with arcpy.da.SearchCursor("TestPolygons", ["SHAPE@", "OID@"]) as s_cursor:
        for shp, oid in s_cursor:
            ic = create_incircle(shp)
            i_cursor.insertRow([ic, oid])

 

JohannesLindner_0-1670429585805.png

 

This is quite slow (1 - 2 seconds per feature), probably because the densification in line 5 adds lots of vertices. This can probably be optimized much, but it works as a quick-and-dirty approach. Sadly, it requires an Advanced license for the Thiessen polygons.


Have a great day!
Johannes

View solution in original post

1 Reply
JohannesLindner
MVP Frequent Contributor

AFAIK, there is no out of the box tool for that.

This thread lists a few algorithms the get an incircle: algorithm - Largest circle inside a non-convex polygon - Stack Overflow

 

I took one of the approaches:

  • tesselate the polygon with a Voronoi diagram (Thiessen polygons in ArcGIS)
  • find the node with the greates distance to the polygon's border, that's the center of the incircle
def create_incircle(polygon):
    add = arcpy.env.addOutputsToMap
    arcpy.env.addOutputsToMap = False
    # densify polygon to get rid of curves, which only have 2 vertices
    polygon = polygon.densify("DISTANCE", 10, 1)
    # extract vertices into a point fc
    vertices = arcpy.management.CreateFeatureclass("memory", "vertices", "POINT", spatial_reference=polygon.spatialReference)
    with arcpy.da.InsertCursor(vertices, ["SHAPE@"]) as cursor:
        for part in polygon:
            for v in part:
                cursor.insertRow([v])
    # create Thiessen polygons
    thiessen = arcpy.analysis.CreateThiessenPolygons(vertices, "memory/thiessen")
    # intersect these with  themselves to get candidate points
    test_points = arcpy.analysis.Intersect([thiessen], "memory/TestPoints", output_type="POINT")
    # get the point inside the polygon with the greatest distance to the polygon's boundary
    candidates = []
    boundary = polygon.boundary()
    with arcpy.da.SearchCursor(test_points, ["SHAPE@"]) as cursor:
        for p, in cursor:
            if p.disjoint(polygon):
                continue
            candidates.append([p, p.distanceTo(boundary)])
    center = sorted(candidates, key=lambda c: c[1])[-1]
    # buffer that point ith the distance to the boundary, return
    arcpy.env.addOutputsToMap = add
    return center[0].buffer(center[1])



# create the ouput fc
incircles = arcpy.management.CreateFeatureclass("memory", "Incircles")
arcpy.management.AddField(incircles, "FID", "LONG")

# create the incircles
with arcpy.da.InsertCursor(incircles, ["SHAPE@", "FID"]) as i_cursor:
    with arcpy.da.SearchCursor("TestPolygons", ["SHAPE@", "OID@"]) as s_cursor:
        for shp, oid in s_cursor:
            ic = create_incircle(shp)
            i_cursor.insertRow([ic, oid])

 

JohannesLindner_0-1670429585805.png

 

This is quite slow (1 - 2 seconds per feature), probably because the densification in line 5 adds lots of vertices. This can probably be optimized much, but it works as a quick-and-dirty approach. Sadly, it requires an Advanced license for the Thiessen polygons.


Have a great day!
Johannes