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).
Good idea. I can see a few use-cases for this. Did you run the second picture as batches and clip?
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. 🙂
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)
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!
@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)
@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)
Well done, both of you. It works like a dream!
Melden Sie sich an, um zu posten, Inhalte zu folgen und mehr. Neu hier? Kostenlos registrieren.
