1

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

Vince
  • 20,017
  • 15
  • 45
  • 64
Wernazuma
  • 75
  • 6

1 Answers1

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:

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.

FelixIP
  • 22,922
  • 3
  • 29
  • 61
  • 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 '20 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 '20 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 '20 at 20:44
  • https://desktop.arcgis.com/en/arcmap/10.3/analyze/creating-tools/adding-a-script-tool.htm – FelixIP Aug 12 '20 at 20:49
  • Thank you so much, I managed to get it going. Good thinking! – Wernazuma Aug 12 '20 at 21:32