Similar to Creating polygons from points using lines as barriers, which has not been resolved; plus, the scenario is slightly different.

I have one dataset with 60,000 points and one with boundaries of ca. 5000 districts. I need to establish areas for the points (using Thiessen). However, the resulting areas must not cross district boundaries.

As a first step, I'd just created the Thiessen polygons. Then I'd make a union with the districts. I could also delete (or dissolve) the residuals that contain no points. (in the example image: black lines)

But: this actually solves nothing. In the example image, the residual of the western district (encircled in blue) would have to be divided up between points A and B. I can't think of any solution to this. example of Thiessen/district combination

1 Answer 1


Script below will work if your points and polygons are stored in shapefiles and have the same defined projection. You also need to run Near tool between points and polygons.

import arcpy
# naming and environment
from arcpy import env
env.workspace = "in_memory"
points,pgons = "points","pgons"
insidePoints, voronois = "pointsInside","voronois"
outFC = "C:/scratch/Thiessens.shp"
# create list of polygons
g = arcpy.Geometry()
geomList = arcpy.management.CopyFeatures(pgons,g)

# handle unique polygons            
def getThisPolygon(i):
    clipper = geomList[i]
    env.extent = clipper.extent
    arcpy.analysis.Select(points, insidePoints, '"NEAR_FID" =%i'%i)
    arcpy.analysis.CreateThiessenPolygons(insidePoints, voronois)
    clipped = "clipped%i" %i
    arcpy.analysis.Clip(voronois, clipper, clipped)
    arcpy.AddMessage("%i completed" %(i+1))
    return clipped

rList = map(getThisPolygon, range(len(geomList)))
env.extent = "MAXOF"
arcpy.management.Merge(rList, outFC)

Original polygons shown by bold black line:

enter image description here

You'll need to change name of output and also transfer point attributes to proximity polygons, unless you'd like to change line where arcpy creates them to include all attributes of the points.

  • I suppose this script is meant to run from inside ArcMap somehow? I tried to define the input and temporary voronois using absolute paths, but somehow I always end up with an empty output shapefile. Already the temporary voronoi shapefile is empty. I only changed some variables: pgons = "D:/DTL/Shapefiles/c27_bounding.shp" points = "D:/DTL/Shapefiles/VID_27.shp" insidePoints = "pointsInside" voronois = "D:/DTL/Shapefiles/temp/voronois.shp" outFC = "D:/DTL/Shapefiles/Thiessen27.shp"
    – Wernazuma
    Aug 12, 2020 at 13:46
  • PS: I did run the near analysis, and I did check that the projection was the same in both files.
    – Wernazuma
    Aug 12, 2020 at 14:11
  • Yes run it from mxd Name your inputs as they are in script. Make sure projection of dataframe is the same as inputs env. Extent = union of inputs or default.
    – FelixIP
    Aug 12, 2020 at 20:44
  • Thank you so much, I managed to get it going. Good thinking!
    – Wernazuma
    Aug 12, 2020 at 21:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.