Calculate the number of fixed-size circles which can fit in an arbitrary polygon (Part I)

[VinceAngelo]

I've seen a number of queries recently, trying to find a tool which would allow a user to determine the maximum number of circles of a certain size which fit inside an arbitrary polygon.  There is, of course, no such COTS tool, but generating a custom tool task starts with generating a region of tightly packed circles, and there's almost a tool to do that, though it needs some tweaking...

The ArcGIS GenerateTessellation tool (arcpy.management.GenerateTessellation) will generate tessellations of a number of given shapes*:

*Note: Transverse Hexagon and Diamond only available in ArcGIS Pro

but "Circle" isn't included because the definition of tessellation ("a covering of an infinite geometric plane without gaps or overlaps by congruent plane figures of one type or a few types" -- Merriam-Webster) requires the figures to be "without gaps or overlaps".

However, if your intent is to use circles anyway, then there are clearly two base geometries one could use to start: SQUARE (which is the same as DIAMOND, just rotated 45 degrees), and HEXAGON (for which TRANSVERSE_HEXAGON is the rotated form, at 90 degrees).  Since the goal here is "packing", the choice should clearly be for HEXAGON (since they are the bestagons )

And how do you go from hexagon to circle? Ahh, well there are a couple of issues to consider:

1. Regular hexagons are composed of six equilateral triangles, so "side" and "radius" are actually the same
2. The radius of a hexagon is the radius of the circumscribing circle, not the inscribed circle, and circumscribed circles of neighboring tiles would overlap (doh!)
3. The Generate Tessellation tool uses hexagon area to determine hexagon size, not side length/radius*
   
   *This is cool, though, because there's a formula to derive area ( A ) from side ( a )
   

So what's the key thing to know? In order to use hexagons to generate packed circles, we need to define the hexagon with respect to the midpoint of the hexagon side (where the inscribed circle touches). 

And what is the difference? I'll save the math, for now, and cut to the chase: You need to increase the circle radius by the inverse of the sine of the interior angle size of a hexagon -- 60 degrees (aka pi / 3 radians).  In Python, that works out to:

import math

diameter = 100
radius   = diameter * 0.5
side     = radius / math.sin(math.pi / 3.0)
hex_area = 1.5 * math.sqrt(3) * side * side

The rest of the problem is simple enough to script:

1. Generate Tessellation with the compensated area
2. Extract centroids from the hexagons as points
3. Use the points to buffer circles of the original diameter

And that would look like this (with an arbitrary Cartesian CRS and extent, and kludged naming):

sr = arcpy.SpatialReference(3968)
tess_shape = 'TRANSVERSE_HEXAGON' # or 'HEXAGON'
arcpy.env.workspace = r"C:\Temp\packed_circles\cartesian.gdb"
extent = '-1000 -1000 1000 1000'
subscript  = 1

print("Generating hexagons with side {:.6f}...".format(side))
arcpy.management.GenerateTessellation(
    Output_Feature_Class="hexagons{:d}".format(subscript),
    Extent=extent,
    Shape_Type=tess_shape,
    Size="{:f} SquareMeters".format(hex_area),
    Spatial_Reference=sr.exportToString()
)
print("Generating centroids...")
arcpy.management.FeatureToPoint(
    in_features="hexagons{:d}".format(subscript),
    out_feature_class="centroids{:d}".format(subscript),
    point_location="CENTROID"
)
print("Generating circles (radius {:.1f})...".format(radius))
arcpy.analysis.Buffer(
    in_features="centroids{:d}".format(subscript),
    out_feature_class="circles{:d}".format(subscript),
    buffer_distance_or_field="{:.1f} Meters".format(radius),
    line_side="FULL",
    line_end_type="ROUND",
    dissolve_option="NONE",
    dissolve_field=None,
    method="PLANAR"
)

Which gives us 100m diameter circles, packed tightly:

And how do I use these circles to find the maximum number of circles in an arbitrary polygon?  Well, that's going to be a bit trickier (aka, "I haven't coded it yet"), so I'll address it in Part II.

- V