Author Archives: dkwiens

Query Planetlabs API from SlackHQ

Lately I’ve been looking into microservices and webhooks, and along the way made a slash command in Slack that will query the Planet image catalog and return matching image IDs. Here’s how:

1.) Make yourself a microservice. I used because it’s free, and I’m cheap. The code I used is here. Basically, it receives a request from Slack (including your api key, item type, asset type, cloud cover limit, and bounding box), formats it into a Planet-friendly query, sends the request to Planet, and returns the response to Slack when it’s ready.

2.) Create a new slash command in Slack. Make sure you set the request URL to your microservice above. I called my slash command ‘/planetcatalog’, but you can use whatever you want.

3.) Install your slash command in Slack.

4.) Run the command, like so (this isn’t a real api key – you’ll need your own):

/planetcatalog a44eb36fe3854f1db1ebXXXdff1a1a3f | PSScene4Band | analytic_sr | 0.1 | {u'type': u'Polygon', u'coordinates': [[[-118.364429, 34.0597], [-118.364429, 34.061552], [-118.362153, 34.061552], [-118.362153, 34.0597], [-118.364429, 34.0597]]]}

The result should be a list of image IDs. Enjoy!


The Trees of Prince George


I’ve been diving into ArcGIS Pro/ArcGIS Online lately, and here’s a quick way to produce a 3D web scene from publicly available data (final example: here). You do need an ArcGIS Online Organizational account with permission to publish in order to follow along.

  1. Find your data. In my case, I used a point feature class representing city trees, provided by the City of Prince George, through their Open Data Catalogue (layer here).
  2. In ArcGIS Pro, start a new local scene project and add the trees as a preset layer (Map -> Add Preset -> Realistic trees. Convenient!). Customize the preset symbology to use GenusName as the type, TreeHeight as the height, proportional as the crown width, and meters as the unit.
  3. Eventually, we need to share the trees as a web feature layer. Unfortunately, web feature layers are limited to 2000 features (there are 3928 trees in the dataset), so we need to split it up. You can do that by making a copy of the feature layer (copy/paste), and setting a definition query for each one (e.g. FID is Less Than 2000/FID is Greater Than or Equal to 2000).
  4. It’s time to share the web scene – Share -> Share As -> Web Scene.
  5. Fill out the web scene metadata (e.g. name, summary, tags, etc.).
  6. Click Analyze to make sure all is well, and troubleshoot the problems as necessary. You may need to change the scene coordinate reference system to match the layers.
  7. Click Share and the web scene should be uploaded to your ArcGIS organizational account. Make sure you share with everyone if this is meant to be a public scene.
  8. In a browser, log into your ArcGIS organizational account. You should see a web scene, hosted feature layer, and service definition. You can customize the web scene if you’d like.
  9. Click on the web scene item. Click Create Web App -> Using a Template. I used the simple scene viewer template. The rest takes you through some configuration steps, but they isn’t a whole lot to configure with this template.
  10. That’s it! You can see my example here:

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!

Basic Mapzen/Tangram test


Over the past few days, I’ve jumped into the Mapzen/Tangram ecosystem. First impressions: somewhat steep learning curve to grasp YAML format, lighting options are pretty neat, very fast to render GeoJSON, documentation is good although asks a lot of the newbie reader, the provided Mapzen Vector Tile Service (mostly OSM) is awesome, and I bet shaders are cool, but I don’t understand them (yet).

My example is largely based off the Mapzen JS Walkthrough, with a few simple additional wrinkles. I draw parcel geometry and attribute data directly from Open Data Prince George. I take building data from the same source, remove OSM buildings from the scene, and draw/extrude 3D buildings based on the building height attribute.