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:
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')

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]
poly_boundary = poly.boundary()
dist = q_dist[2] if q_dist[3] else -q_dist[2]
max_dist = 0
dists = {}
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%
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)
for i in range(10): # quadify 10x
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.

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!