The main script associated with a tool. Download the tooset from Triangulation Tools for use in ArcMap versions 10.0 and above
'''
Vor_Del_Main.py
Author
Dan Patterson
Dept of Geography and Environmental Studies
Carleton University, Ottawa, Canada
Dan_Patterson@carleton.ca
Original Delaunay/Voronoi code written by Steven Fortune
Translation to Python by Bill Simons
A variant also exists which uses the same Delaunay code,
http://arcscripts.esri.com/details.asp?dbid=16684
by
David Bissessar
School of Computer Science
Carleton University, Ottawa, Canada
dbissess@connect.carleton.ca
See header comments in vor_del.py for more on the origins and the
original authors websites.
Purpose
Produce Delaunay Triangulation and Voronoi Diagram for point data.
Requires
vor_del.py, arcpy_helper.py installed in the same folder as this and toolbox.
Parameter list
Parameter Properties
Display Name Data type Type Direction MultiValue Dependency
argv[1] Point layer Feature layer Required Input No
argv[2] Output Shapefile Feature class Required Output No
argv[3] Output type String Required Input No
Filter Value list consisting of Polygon, Polyline
Default Polygon
argv[4] Form String Required Input No
Filter Value list consisting of Delaunay, Voronoi
Default Delaunay
'''
#-----------------------------------------------------------------------------------------------------------------------------------------------------------
degToRad = math.pi/180.0
def circleMake(aPnt,aRadius):
'''create a circle from a point [X,Y] and radius
returns a list of points representing the circle as an Ngon
'''
X0, Y0 = aPnt[0], aPnt[1]
X1, Y1 = X0 + aRadius, Y0
circPnts = [[X1, Y1]]
for i in range(1,360):
X = X0 + math.cos(degToRad * i) * aRadius
Y = Y1 + math.sin(degToRad * i) * aRadius
circPnt = [X,Y]
circPnts.append(circPnt)
del circPnt
return circPnts
def extentCenter(pnts):
'''Returns the average and median X, Y coordinates of a series
of input points'''
N = len(pnts)
Xsum = 0.0; Ysum = 0.0
ReturnList = []; xList = []; yList = []
for pnt in pnts:
X = pnt[0]; Y = pnt[1]
xList.append(pnt[0]); yList.append(pnt[1])
L = min(xList); R = max(xList); B = min(yList); T = max(yList)
Xcent = L + (R - L)/2.0; Ycent = B + (T - B)/2.0
return [Xcent, Ycent]
def extentPnts(pnts):
'''Returns the min, max X, Y of input points'''
N = len(pnts)
Xsum = 0.0; Ysum = 0.0
ReturnList = []; xList = []; yList = []
for pnt in pnts:
# = pnt[0]; Y = pnt[1]
xList.append(pnt[0]); yList.append(pnt[1])
L = min(xList); R = max(xList); B = min(yList); T = max(yList)
#print "N ", str(N), "Xcent", str(Xcent),"Ycent ", str(Ycent)
LL = [L,B]; UL = [L,T]; UR = [R,T]; LR = [R,B]
return [LL, UL, UR, LR]
def buffer_extent(extent_list, buffer_dist):
'''produce a square buffer from the LL, UL, UR, LR points'''
LL, UL, UR, LR = extent_list
dx = buffer_dist
dy = buffer_dist
L = LL[0] - dx; R = LR[0] + dx
T = UL[1] + dy; B = LL[1] - dy
LL = [L,B]; UL = [L,T]; UR = [R,T]; LR = [R,B]
return [LL, UL, UR, LR]
def infinity_circle(pnts):
'''create an infinity circle to append to the original point list'''
cent_pnt = extentCenter(pnts)
LL, UL, UR, LR = extentPnts(pnts)
dx = LR[0] - LL[0]
dy = UL[1] - LL[1]
circle_radius = max([dx, dy]) * 10 #increase radius by a factor of 10
circPnts = circleMake(cent_pnt, circle_radius)
return circPnts
def list_to_pnts(pnts):
'''convert a list to arcpy point objects'''
out_pnts = []
for i in pnts:
out_pnts.append(arcpy.Point(i[0],i[1]))
return out_pnts
def points_to_list(pnts):
'''convert point objects to a list of lists [[1,1],[2,2]]'''
aList = []
for pnt in pnts:
aList.append([pnt.X, pnt.Y])
return aList
def pnts2polyline(vertices, indices):
'''create a line from 2 points, format: edge index, pnt1 index, pnt2 index'''
merged = []
for triplet in indices:
if -1 not in triplet:
vertexPoints = []
start_pnt = list(vertices[triplet[1]])
end_pnt = list(vertices[triplet[2]])
vertexPoints.append(start_pnt)
vertexPoints.append(end_pnt)
merged.append( [vertexPoints ] ) #need the extra list
return merged
def triangle2polyline(vertices, indices):
'''create a closed-loop polyline from 3 points'''
merged = []
for trio in indices:
vertexPoints = []
for vertexIndex in trio:
x = vertices[vertexIndex].X
y = vertices[vertexIndex].Y
vertexPoints.append( [x,y] )
vertexPoints.append(vertexPoints[0])
merged.append( [vertexPoints ] ) #need the extra list
return merged
def triangle2polygon(vertices, indices):
'''create a polygon from 3 points'''
merged = []
for trio in indices:
vertexPoints = []
for vertexIndex in trio:
x = vertices[vertexIndex].X
y = vertices[vertexIndex].Y
vertexPoints.append( [x,y] )
merged.append( [vertexPoints ] ) #need the extra list
return merged
def sort_radial(pnts, cent):
'''performs a radial sort about the center point'''
sorted_pnts =[]
out_pnts = []
for pnt in pnts:
p_angle = (pnt - cent).angle()
sorted_pnts.append([p_angle, pnt])
sorted_pnts.sort()
sorted_pnts.reverse()
for i in sorted_pnts:
out_pnts.append(i[1])
return [out_pnts, sorted_pnts]
#-----------------------------------------------------------------------------------------------------------------------------------------------------------
#Main section
#Import the modules, create the geoprocessor and describe objects
# and get some base information about the feature layer
# fc is the feature class on disk used for spatial ref template
import sys, os, math
import arcpy
import arcpy_helper
#reload(arcpy_helper) #uncomment if editing arcpy_helper
import vor_del
#reload(vor_del) #uncomment if editing vor_del
#arcpy, gp_version = arcpy_helper.gp_create(10)
gp_version = 10
arcpy.AddMessage("\n Voronoi / Delaunay tool. Processing: using version... " + \
str(gp_version))
toolboxes = ["Data Management Tools.tbx", "Analysis Tools.tbx"]
arcpy, msg, passed = arcpy_helper.get_toolboxes(toolboxes, arcpy)
arcpy.AddMessage(msg)
if not passed:
arcpy.AddMessage("\n Exiting ..... \n")
del arcpy
sys.exit()
arcpy.env.overWriteOutput = True #overwrite temporary files
# 1 Obtain the arguments, determine the spatial reference, shape
# type, shapefield and FID field
in_FC = sys.argv[1]
out_FC = sys.argv[2]
output_type = sys.argv[3]
output_form = sys.argv[4]
buffer_dist = float(sys.argv[5])
desc = arcpy.Describe
fc = desc(in_FC).CatalogPath.replace("\\","/")
SR = desc(in_FC).spatialReference
theType = desc(in_FC).ShapeType
shapeField = desc(in_FC).ShapeFieldName
OIDField = desc(in_FC).OIDFieldName
# 2 Read the data and collect into a point list
rows = arcpy.SearchCursor(in_FC)
input_pnts =[]
for row in rows:
a_shape = row.Shape
pnt = a_shape.getPart()
input_pnts.append(pnt)
del rows
del row
# 3 Run the Delaunay or the Voronoi and produce the outputs
# In the case of Voronoi, create an infinity circle so that the tesselation
# is unbounded, then clip to the desired extent if polygon output is chosen
if output_form == "Delaunay":
triangleIndices = vor_del.computeDelaunayTriangulation(input_pnts)
if output_type == "Polyline":
triangleVertices = triangle2polyline( input_pnts, triangleIndices)
else:
triangleVertices = triangle2polygon( input_pnts, triangleIndices)
triangleVertices.sort() #sort left to right, bottom to top
theFields = [] #optional fields to add
arcpy_helper.polyFromPoints(out_FC, output_type, SR, triangleVertices, theFields, arcpy)
else: #Voronoi output
orig_list = points_to_list(input_pnts) #determine the extent of the original input points
extent_pnts = extentPnts(orig_list)
if buffer_dist > 0.0:
extent_pnts = buffer_extent(extent_pnts, buffer_dist)
infinity_pnts = infinity_circle(orig_list)
temp_pnts = list_to_pnts(infinity_pnts)
input_pnts.extend(temp_pnts) #add the infinity circle to the point list
theFields = []
extent_FC = out_FC.replace(".shp", "_extent.shp") #extent feature class
temp_FC = out_FC.replace(".shp", "_temp.shp") #temporary FC
#
result = vor_del.computeVoronoiDiagram(input_pnts)
vert_out = result[0]
#lines_out = result[1] #not used
edges_out = result[2]
triangles_out = result[3]
if output_type == "Polyline": #polyline output
segments = pnts2polyline(vert_out, edges_out)
arcpy_helper.polyFromPoints(extent_FC, "Polygon", SR, [[extent_pnts]], theFields, arcpy)
arcpy_helper.polyFromPoints(temp_FC, output_type, SR, segments, theFields, arcpy)
arcpy.Clip_analysis(temp_FC, extent_FC, out_FC)
try:
arcpy.CalculateField_management(out_FC, "ID", "!FID!", "PYTHON")
except:
pass
else: #polygon output
vert_list, site_dict, poly_list = vor_del.voronoi_polygons(input_pnts, vert_out, triangles_out)
poly_verts = []
theFields = []
for i in poly_list:
site_pnt = i[0]
poly_pnts = i[1]
if vor_del.point_in_poly(site_pnt, poly_pnts):
p_list = points_to_list(i[1])
poly_verts.append([p_list])
#
union_FC = out_FC.replace(".shp","_union.shp") #unioned FC
clip_FC = out_FC.replace(".shp", "_clip.shp") #clipped FC
arcpy_helper.polyFromPoints(extent_FC, output_type, SR, [[extent_pnts]], theFields, arcpy)
arcpy_helper.polyFromPoints(temp_FC, output_type, SR, poly_verts, theFields, arcpy)
arcpy.Union_analysis([extent_FC, temp_FC], union_FC, "NO_FID")
arcpy.Clip_analysis(union_FC, extent_FC, clip_FC)
arcpy.MultipartToSinglepart_management(clip_FC, out_FC) #final output file
del_fields = ["ID", "ID_1"]
try: #try deleting the extraneous fields
arcpy.DeleteField_management(out_FC, del_fields)
except:
pass
del_list = [temp_FC, union_FC, clip_FC]
for i in del_list:
arcpy.AddMessage("Deleting temporary file: " + str(i) )
arcpy.Delete_management(i)
#
#optional point file
temp_FC = os.path.dirname(out_FC).replace("\\","/") + "/vertices.shp"
vert_FC = out_FC.replace(".shp", "_vertices.shp") #temporary FC
arcpy_helper.createPointFile(temp_FC, "Point", SR, [vert_out], theFields, arcpy)
if output_type == "Polygon":
arcpy.Intersect_analysis([temp_FC, extent_FC], vert_FC, "ALL")
arcpy.Delete_management(temp_FC)
arcpy.AddMessage("\nExtra files have also been created showing the extent polygon and the voronoi vertices")
arcpy.AddMessage(" extent: " + str(extent_FC))
arcpy.AddMessage(" vertices: " + str(vert_FC))
arcpy.AddMessage("Processing complete" + "\n")
del arcpy
No comments:
Post a Comment