Monthly Archives: August 2016

Mapbox Visual Center Algorithm: arcpy (attempt #2)


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:
        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,'#')
        label_points.append(arcpy.PointGeometry(arcpy.Point(best_point[0], best_point[1]),sr))
        for cell in best_point[2]:
arcpy.CopyFeatures_management(label_points, r'in_memory\label_points')
arcpy.CopyFeatures_management(cells, r'in_memory\cells')

Mapbox Visual Center Algorithm: arcpy


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
    for k,v in dists.iteritems():
        if k > max_dist * 0.90: # precision = 90%
    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


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!

Simple CesiumJS Skydive Sim


Note: updated version with guidance lights here.

On July 30, 2016, a man named Luke Aikins went skydiving, without a parachute, and landed safely in a net. I was curious what that would feel like, so I made this simple skydiving simulation using CesiumJS.

You can track using W, A, S, D to feel how much control you could possibly have on the descent.

The simulation parameters are as follows:

  • the location is similar to the event, but I’m not exactly sure where it took place. Apparently, it was at the Big Sky movie ranch, near Simi Valley, CA. This is the nearest open area to where Google Maps tells me that is.
  • the red target on the ground is 30m x 30m, same as at the event
  • the simulation starts at 7620m (or 25,000 feet), same as at the event
  • you fall at 53.6448 meters/second, which is 120mph, which is at the slower end of the velocity range a human would travel during free fall. I figure he was trying to go as slow as possible, but there’s only so much you can do to make that happen!
  • I was surprised to learn (simultaneously from both @tkw954 and @erikfriesen – thanks!) that a good skydiver can track at 45 degrees (1:1 glide ratio), which made the math for horizontal movement super-simple – it equals the fall rate.

I stole borrowed most of the animation/control code from the Cesium Camera Tutorial.

The simulation ends in blackness, which is because I was too lazy to make a net, however with no net, that seems fitting. Have fun!

edit (08/06/16): this post was featured on Maps Mania!

edit (08/09/16): added guidance lights to this example. This is likely not how the actual light guidance system worked. In my example, the farther you’re aimed from each corner (up to 100m), the more red the circle turns. The closer you’re aimed to each corner, the more green it turns. It would be more useful to code the lights to indicate in which direction to correct, but I only have so much time.

edit (09/01/16): this post was mentioned in the Cesium Blog!

edit (12/29/16): this post was mentioned in the Exploratorium Blog!