Yesterday, I attempted to replicate the Mapbox Visual Center algorithm using ArcGIS arcpy (Python). Today, I actually looked at the source code (code, license) and attempted to fully port it over from JavaScript. I think it works properly, but you tell me.
def polylabel(polygon, centroid, precision, debug):
cells = []
precision = precision or 1.0
extent = polygon.extent
minX, minY, maxX, maxY = (extent.XMin, extent.YMin, extent.XMax, extent.YMax)
width = extent.width
height = extent.height
cellSize = min([width, height])
h = cellSize / 2
cellQueue = []
x = minX
while x < maxX:
y = minY
while y < maxY: cell = Cell(x+h, y+h, h, polygon) cells.append(cell.geom) cellQueue.append(cell) y += cellSize x += cellSize bestCell = Cell(centroid[0], centroid[1], 0, polygon) numProbes = len(cellQueue) while len(cellQueue): cellQueue.sort(key=lambda cell: cell.d) cell = cellQueue.pop() cells.append(cell.geom) if cell.d > bestCell.d:
bestCell = cell
if cell.max - bestCell.d <= precision:
continue
h = cell.h/2
cellQueue.append(Cell(cell.x-h, cell.y-h, h, polygon))
cellQueue.append(Cell(cell.x+h, cell.y-h, h, polygon))
cellQueue.append(Cell(cell.x-h, cell.y+h, h, polygon))
cellQueue.append(Cell(cell.x+h, cell.y+h, h, polygon))
numProbes += 4
return [bestCell.x, bestCell.y, cells]
class Cell():
def __init__(self, x, y, h, polygon):
self.x = x
self.y = y
self.h = h
self.d = pointToPolygonDist(x, y, polygon)
self.max = self.d + self.h * math.sqrt(2)
self.geom = arcpy.Polygon(arcpy.Array([arcpy.Point(x-h,y-h),arcpy.Point(x+h,y-h),arcpy.Point(x+h,y+h),arcpy.Point(x-h,y+h)]))
def pointToPolygonDist(x, y, polygon):
point_geom = arcpy.PointGeometry(arcpy.Point(x, y),sr)
polygon_outline = polygon.boundary()
min_dist = polygon_outline.queryPointAndDistance(point_geom)
sign_dist = min_dist[2] * ((min_dist[3]-0.5)*2)
return sign_dist
fc = 'BEC_POLY selection'
sr = arcpy.Describe(fc).spatialReference
centroids = []
label_points = []
cells = []
with arcpy.da.SearchCursor(fc,['SHAPE@','SHAPE@XY']) as cursor:
for row in cursor:
best_point = polylabel(row[0],row[1],.75,'#')
centroids.append(arcpy.PointGeometry(row[0].centroid,sr))
label_points.append(arcpy.PointGeometry(arcpy.Point(best_point[0], best_point[1]),sr))
for cell in best_point[2]:
cells.append(cell)
arcpy.CopyFeatures_management(centroids,r'in_memory\centroids')
arcpy.CopyFeatures_management(label_points, r'in_memory\label_points')
arcpy.CopyFeatures_management(cells, r'in_memory\cells')
Advertisements
