Voronoi partitions (Thiessen Polygons) with limited extents from a multi polygon dataset.

[CliveCartwright]

I feel that the Voronoi partitions (Thiessen Polygons) created using the Create Thiessen Polygons tool would benefit with an added functionality that allows the user to restrict the calculation of the Thiessen Polygons to the boundaries of another dataset (e.g. district boundaries). 

[DavidPike]

Good idea.  I can see a few use-cases for this.  Did you run the second picture as batches and clip?

[CliveCartwright]

Hi @DavidPike, 

Sadly no. I manually mocked it up as helpful visual to get the idea across. But now you've said that, I'm tempted to test my mock by doing just as you've suggested. 🙂 

[BrennanSmith1]

This seems like a useful feature, so I scripted a workaround.

import arcpy
arcpy.env.addOutputsToMap = False
arcpy.env.overwriteOutput = True

# Input feature classes for points and polygon zones. Also the output
boundary_fc = "v_poly"
points_fc   = "v_point"
output_fc   = "v_output"


thiessen_feature_classes = []
# loop through poylgon OIDs
with arcpy.da.SearchCursor(boundary_fc, ["OID@"]) as boundary_cursor:
    for row in boundary_cursor:
        oid = row[0]

        # Create feature layers
        arcpy.management.MakeFeatureLayer(points_fc, "points_lyr")
        arcpy.management.MakeFeatureLayer(boundary_fc, "boundary_lyr")

        # Select the current boundary polygon
        arcpy.SelectLayerByAttribute_management("boundary_lyr", "NEW_SELECTION", f"OBJECTID = {oid}")

        # Select points within the selected boundary polygon
        temp_points_fc = r"memory\temp_points"
        arcpy.SelectLayerByLocation_management("points_lyr", "WITHIN", "boundary_lyr", "", "NEW_SELECTION")
        arcpy.management.CopyFeatures("points_lyr", temp_points_fc)

        # Check if any points were selected
        if int(arcpy.GetCount_management(temp_points_fc).getOutput(0)) > 0:
            # Create intermediate thiessen polygons
            output_thiessen_fc = r"memory\Thiessen_OID_" + str(oid)
            arcpy.env.extent = boundary_fc
            arcpy.analysis.CreateThiessenPolygons(temp_points_fc, output_thiessen_fc)

            # Clip the Thiessen polygons to the boundary
            clipped_thiessen_fc = r"memory\Clipped_Thiessen_OID_" + str(oid) 
            arcpy.analysis.Clip(output_thiessen_fc, "boundary_lyr", clipped_thiessen_fc)

            # Keep a running list of these clipped thiessen polys
            thiessen_feature_classes.append(clipped_thiessen_fc) 
            
        #remove intermediate layers
        arcpy.management.Delete("points_lyr")
        arcpy.management.Delete("boundary_lyr")

# merge them all
if thiessen_feature_classes:
    arcpy.env.addOutputsToMap = True
    arcpy.management.Merge(thiessen_feature_classes, output_fc)
    

[DavidPike]

Very nice.  I wonder if it can be faster by taking the arcpy.management.MakeFeatureLayer() out of the cursor for points and boundaries.  they might be instantiated in each loop unnecessarily, but could well be misinterpreting!

[BrennanSmith1]

@DavidPike good idea! I got frustrated when I couldn't get a scipy.spatial.Voronoi approach to work, so I just slapped it together using the GP tools. Here is a slightly optimized approach:

import arcpy
arcpy.env.addOutputsToMap = False
arcpy.env.overwriteOutput = True

# Input feature classes for points and polygon zones, and the output polygon. Assuming we're working in a gdb workspace.
boundary_fc = "v_poly"
points_fc   = "v_point"
output_fc   = "v_output"

# create feature layers for selections
points_lyr   = arcpy.management.MakeFeatureLayer(points_fc, "points_lyr")
boundary_lyr = arcpy.management.MakeFeatureLayer(boundary_fc, "boundary_lyr")

# get polygon oids and extent
oid_field = arcpy.Describe(boundary_fc).OIDFieldName
oids = [row[0] for row in arcpy.da.SearchCursor(boundary_fc, ["OID@"])]
arcpy.env.extent = boundary_fc

# empty list of thiessen FCs to merge at the end
thiessen_feature_classes = []

# loop through poylgon OIDs
for oid in oids:
    
    # Select the current boundary polygon
    arcpy.management.SelectLayerByAttribute(boundary_lyr, "NEW_SELECTION", f"{oid_field} = {oid}")

    # Select points within the selected boundary polygon
    arcpy.management.SelectLayerByLocation(points_lyr, "WITHIN", boundary_lyr, "", "NEW_SELECTION")
    
    # Check if any points were selected
    if arcpy.Describe(points_lyr).FIDSet:
        
        # Create intermediate thiessen polygons
        temp_fc = r"memory\Thiessen_OID_" + str(oid)
        arcpy.analysis.CreateThiessenPolygons(points_lyr, temp_fc)

        # Clip the Thiessen polygons to the boundary
        clipped_fc = r"memory\Clipped_Thiessen_OID_" + str(oid) 
        arcpy.analysis.Clip(temp_fc, boundary_lyr, clipped_fc)

        # Keep a running list of these clipped thiessen polys
        thiessen_feature_classes.append(clipped_fc) 
        
#remove intermediate layers
arcpy.management.Delete("points_lyr")
arcpy.management.Delete("boundary_lyr")

# merge them all
if thiessen_feature_classes:
    arcpy.env.addOutputsToMap = True
    arcpy.management.Merge(thiessen_feature_classes, output_fc)

[DavidPike]

@BrennanSmith1 that looks perfect.  Just my final recommendations for posterity as I'm sure people will copy your script for future use:

# Check if any points were selected
# if arcpy.Describe(points_lyr).FIDSet:
### selectlayerByLocation has a derived output for 
### number of features selected already.  The describe
### operation is quite expensive I think so perhaps this:
layer, count = arcpy.management.SelectLayerByLocation(points_lyr, "WITHIN", boundary_lyr, "", "NEW_SELECTION")

if int(count)>0:

### may sometimes be issues if boundaries have greater extent
### than the points, resulting in a rectangular clip before reaching
### the boundaries.  But this is thinking from ArcMap days, maybe a long
### distant memory in Pro.    

# Set the extent environment using a feature class 
#(can also maybe be set outside loop??)
arcpy.env.extent = boundary_fc   
# Create intermediate thiessen polygons
temp_fc = r"memory\Thiessen_OID_" + str(oid)
arcpy.analysis.CreateThiessenPolygons(points_lyr, temp_fc)

[CliveCartwright]

Well done, both of you. It works like a dream!