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.

Antipodes (MapboxGL)

In 2013, I made an interactive antipode map (points on opposite sides of the Earth) using the Google Maps API. As I’m trying to make the leap from Google Maps/Mapbox JS to Mapbox GL JShere is my rehash of the same map. The only tangible difference between these two maps is that you can rotate this version (right-click and drag).

Some issues I have yet to solve:

• I can’t seem to get the rotated maps’ labels to display correctly (i.e. right-side up and left to right). I would need to somehow change them to be rendered as a mirror image.
• The second set of maps don’t render at all in Chrome, although they do while testing locally. This entire project seems to work in IE and Firefox. edit: this project now seems to work in Chrome.

Draggable Latitude Line Map

I just returned home from a trip to Scandinavia (Denmark/Sweden/Norway), and was curious to see how far north I’d gone. I made it as far as Tromsø, which is at about 69.7°N, but where would that latitude place me in North America? It turns out, pretty far north – about halfway up Baffin Island. See for yourself on this quick map (drag the marker to move the line north and south, then see where in the world passes through that latitude).

JavaScript: UTM ↔ LatLng

I’m sure this has been done (more thoroughly) many times before, but I wanted to take a crack at the simplified UTM/LatLng conversion equations (as listed here). You can find my functions, latLngToUTM(lat,lng) and utmToLatLng(easting,northing,zone), here.

Old Prince George Floodplain Map

Just thought I’d share a cool looking old map overlay I found digging through the BC Map Services. I enjoy the very detailed, yet simple look. I’ve also superimposed the City of Prince George Flood Plain 1997 layer, which I think is the digitized version of the scanned BC map. Kinda scary to see how far up Patricia/Winnipeg/17th/Massey it would flood, but I guess the City’s got it figured out.

WMS GetFeatureInfo Request

Continuing on with yesterday’s exploration into GeoMet and its myriad web mapping services (WMS), today I added GetFeatureInfo requests (map here).

If you’ve used WMS, you’re probably familiar with GetCapabilities (retrieves metadata) and GetMap (retrieves map tiles) requests. GetFeatureInfo is another type of request which retrieves attribute information at a given location, provided the server (e.g. GeoMet) has been configured to provide this information (look for “queryable=’1′ ” in the layer metadata).

One stumbling block I encountered (and likely handled poorly) was overcoming the issue of cross-origin resource sharing (CORS). No matter what I tried, I was blocked from retrieving the request response because, of course, my webpage resides on a different server than the WMS. I got around this by calling a PHP script that copies the text response from the request locally and reads that version. Honestly, I’m not sure why that works and reading the original request doesn’t, but it does and that’s good enough for me. If you have a simpler solution (I’m sure it’s out there), I’d love to learn about it!