Hard line-barriers for Thiessen polygons

Question

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.

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)))
arcpy.AddMessage("Merging...")
env.extent = "MAXOF"
arcpy.management.Merge(rList, outFC)
arcpy.management.Delete("in_memory")

Original polygons shown by bold black line:

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.