Tag Archives: algorithm

Mapbox Visual Center Algorithm: arcpy

quadtree

On Monday, Mapbox published a JavaScript implementation of a fast algorithm for finding the visual center of a polygon using quadtrees (blog, code). I took a stab at it using Python/arcpy in ArcGIS:

def quadify(poly): # split polygon into quads
    ext = poly.extent
    TM = arcpy.Point((ext.XMin+ext.XMax)/2,ext.YMax)
    LM = arcpy.Point(ext.XMin,(ext.YMin+ext.YMax)/2)
    BM = arcpy.Point((ext.XMin+ext.XMax)/2,ext.YMin)
    RM = arcpy.Point(ext.XMax,(ext.YMin+ext.YMax)/2)
    TL = arcpy.Polygon(arcpy.Array([ext.upperLeft,LM,poly.centroid,TM]),sr)
    TR = arcpy.Polygon(arcpy.Array([TM,poly.centroid,RM,ext.upperRight]),sr)
    BL = arcpy.Polygon(arcpy.Array([LM,ext.lowerLeft,BM,poly.centroid]),sr)
    BR = arcpy.Polygon(arcpy.Array([poly.centroid,BM,ext.lowerRight,RM]),sr)
    return [TL,TR,BL,BR]
def inspect(quad,poly): # calculate quad radius & quad center to polygon distance
    poly_boundary = poly.boundary()
    quad_center = arcpy.PointGeometry(quad.centroid,sr)
    radius = quad_center.distanceTo(quad.firstPoint)
    q_dist = poly_boundary.queryPointAndDistance(quad_center)
    dist = q_dist[2] if q_dist[3] else -q_dist[2]
    return dist + radius
def loop_quads(quads): # evaluate the quads & return if in the top 10% of values
    max_quads = []
    max_dist = 0
    dists = {}
    for quad in quads:
        cur_dist = inspect(quad,row[0])
        dists[cur_dist] = quad
        max_dist = cur_dist if cur_dist > max_dist else max_dist
        out_polys.append(quad)
    for k,v in dists.iteritems():
        if k > max_dist * 0.90: # precision = 90%
            max_quads.append(v)
    return max_quads
fc = 'wetlands_select' # feature class/layer
sr = arcpy.Describe(fc).spatialReference # spatial reference
out_polys = []
out_points = []
with arcpy.da.SearchCursor(fc,'SHAPE@',spatial_reference=sr) as cursor: # loop polygons
    for row in cursor:
        ext = row[0].extent
        center = arcpy.PointGeometry(arcpy.Point((ext.XMin+ext.XMax)/2,(ext.YMin+ext.YMax)/2),sr)
        radius = max(center.distanceTo(ext.upperLeft),center.distanceTo(ext.upperRight),center.distanceTo(ext.lowerLeft),center.distanceTo(ext.lowerRight))
        circle = center.buffer(radius)
        quads = quadify(circle) # start with global circle
        for i in range(10): # quadify 10x
            max_quads = loop_quads(quads) # evaluate current set of quads
            quads = []
            for quad in max_quads:
                quads += quadify(quad)
        out_points.append(arcpy.PointGeometry(quad.centroid,sr)) # return one of the best quad centers
arcpy.CopyFeatures_management(out_polys,r'in_memory\polys')
arcpy.CopyFeatures_management(out_points,r'in_memory\points')

Notes:

I’m not 100% sure that this emulates the base Mapbox algorithm, and I didn’t attempt to implement the “priority queue” enhancement.

The results are somewhat sensitive to the precision factor, but I suppose that’s the be expected with a speed algorithm.

I realize that you can calculate a label point much easier within ArcGIS using out-of-the-box tools. This was just me exploring the algorithm.

Anyhow, I welcome your comments on the code above!

Advertisements