perrygeo

About | Articles | CV

Creating UTFGrids directly from a polygon datasource

Mon 20 August 2012

We've begun to rely on the interactivity provided by UTFGrids in many of our recent web maps. (Quick recap: UTFGrids are "invisible" map tiles that allow direct interactivity with feature attributes without querying the server.) Earlier this year, I created the initial OpenLayers UTFGrid support and was glad to see it accepted into OpenLayer 2.12 (with some enhancements).

With the client-side javascript support in place, the only missing piece in the workflow was to create the UTFGrid .json files. We had expirimented with several alternate UTFGrid renderers but Mapnik's rendering was by far the fastest and produced the best results. Using Tilemill was a convenient way to leverage the Mapnik UTFGrid renderer but it came at the cost of a somewhat circuitious and manual workflow:

Load the data up into Tilemill,
Configure interactivity fields
Export to .mbtiles
Convert to .json files

What we really needed was a script to take a polygon shapefile and render the UTFGrids directly to files. Mapnik would provide the rendering while the Global Map Tiles python module would provide the logic for going back and forth between geographic coordinates and tile grid coordinates. From there it's just a matter of determining the extent of the data set and, for a specified set of zoom levels, looping through and using Mapnik to render the UTFGrid to a .json file in Z/X/Y.json directory structure.

Get `create-utfgrids` on github

If we have a mercator polygon shapefile of ecoregions and want to render UTFGrids for zoom levels 3 through 5 using the dom_desc and div_desc attributes, we could use a command like

$ ./create_utfgrids.py test_data/bailey_merc.shp 3 5 ecoregions -f dom_desc,div_desc

WARNING:
This script assumes a polygon shapefile in spherical mercator projection.
If any of these assumptions are not true, don't count on the results!
 * Processing Zoom Level 3
 * Processing Zoom Level 4
 * Processing Zoom Level 5

and inspect the output (e.g. zoom level 5, X=20, Y=18)

$ cat ecoregions/5/20/18.json | python -mjson.tool
{
    "data": {
        "192": {
            "div_desc": "RAINFOREST REGIME MOUNTAINS", 
            "dom_desc": "HUMID TROPICAL DOMAIN"
        }, 
...
    "grid": [
        "  !!!!!!!!!#####$%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%", 
...

Some caveats:

This currently only works for polygon datasets in a Web Mercator projection.
It's only tested with shapefiles as it assumes a single-layer datasource at the moment. Full OGR Datasource support would not be too difficult to add for PostGIS, etc.
It assumes a top-origin tile scheme (as do OSM and Google Maps). Supporting TMS bottom-origin schemes in the future should be straightforward.
Requires OGR and Mapnik >= 2.0 with python bindings. Finding windows binaries for the required version of Mapnik may be difficult so using OSX/Linux is recommended at this time.

Many thanks to Dane Springmeyer for his help on UTFGrid related matters and and to Klokan Petr Přidal for his MapTiler docs

Introducing the Madrona framework

Wed 11 July 2012

Madrona: A software framework for effective place-based decision making

My work at Ecotrust mainly revolves around creating web-based spatial analysis tools - software to bring data-driven science to the place-based descision making process. This began several years ago when I joined the MarineMap team. Since working with Ecotrust, we've taken the MarineMap software far beyond it's original niche. What was once a specific tool for marine protected area planning has now become a powerful framework for all sorts of web-based spatial tools in the realms of marine, forestry, conservation planning, aquatic habitat restoration, etc. So, in a sense, Madrona is a recognition of that evolution.

From the official Madrona release announcement from the Ecotrust blog post:

Over the last year we’ve distilled the best ideas from our most successful tools into a suite of software building blocks that can be mixed and matched to create cutting-edge software for decision support and spatial planning at any scale. These building blocks are already at the heart of our work and now we’re ready to share them with you.

So what is Madrona from a developer's perspective?

A set of python django apps that provide models, views and templates for representing spatial features and solving problems specific to spatial decision tools.
A RESTful API for accessing spatial features
A collection of javascript libraries (based on JQuery) to provide a web-based interface to the API.

In short, we think its a great platform for spatial tools and we want to open it up to the wider developer audience. Ecotrust already has many madrona-based apps in the wild (with many more in development) but we're hoping to get other folks using (and contributing to) the Madrona framework in the future.

I know this post is short on technical details but there will more to come ... for now, check out the technology page for an overview or the developer's page to dive in.

Migrating from Wordpress to Jekyll

Sat 28 April 2012

I just switched this blog from an ancient version of wordpress running on a VPS to a static-file jekyll bootstrap site (hosted by github). Let me know if you experience any wierdness on the site or feeds. I've taken good measures to make sure links don't break (old URLS should get a 301 permanent redirect to blog.perrygeo.net) but let me know if you get any 404s.

So why do it?

Having a PHP-MySQL app running on a VPS just to serve up a bunch of blog posts seemed excessive. I don't have the desire to maintain that sort of infrastructure for a simple blog!
Wordpress' editing and admin interface suck. I prefer vim and bash.
Markdown is a great language for quickly banging out blog posts.
Static files just make sense for what is basically static content.
Github pages provides the hosting for me and even handles CNAMEs for DNS.
Managing revisions with git.

The conversion process

It was not an entirely smooth transition, most of which can be traced directly to dumb decisions on my part. I won't recount the entire process (there are plenty of guides on internets) but I'll outline the major steps here:

Export the wordpress blog to an xml file. I has to use xmllint to clean it up a bit.
Set up a disqus account and import my wordpress file. Disqus will handle all the comments which are the only dynamic content on the page.
Use exitwp.py to convert the xml to jekyll markdown files. This worked OK. Not great. Tags and formatting did not come through as expected and I had to wrestle the script a bit. Tables were destroyed and some iframes (youtube links) were lost.
Forked Jekyll Bootstrap and brought in my posts.
Started tweaking of css and markdown to get formatting right. Still have a ways to go on this front - let me know if there is any content you'd like me to restore faster than others.
Had to write a little web service to redirect posts; the old blog stupidly used the default wordpress URLS like /wordpress/?p=4 which needed to go to /2010/01/01/blah
My images were all over the place; some I had in wordpress uploads, others on various servers, some were absolute links, others relative. Gathering them all in one place and using some sed-fu to get the paths right was essential.
Retagged some posts - still working on tags.
Set up Google Analytics to track usage.

I think that's about it. There are still some big formatting problems on older posts (mostly due to the fact that I used blockquotes for code). And tables are still destroyed. I'll be working on cleaning these up as I go along.

Overall impression of Jekyll-Bootstrap and hosting with Github pages? Awesome. I would highly recomend it to anyone starting a new blog or converting a smaller/better-behaved wordpress site. It is so much better than having to deal with PHP and MySQL (hopefully the last time I'll ever see them!). But the conversion was a bit tricky and took way more of my Friday and Saturday than I'd like to admit. I would not want to do that again... But I'm glad did.

What do you think of the new digs?

Working with mbtiles in python

Sun 25 March 2012

python-mbtiles. Check it out.

I've been working a bit with Tilemill lately and love the Carto css styling, iteractivity through UTFGrids and being able to export the whole deal as a single mbtiles sqlite database. But when it comes to working with the mbtiles databases, I've found both Tilestache and Tilestream to be fairly limiting:

Tilestache serves images but does not (yet) serve up UTFGrids _directly from mbtiles _ while Tilestream hardcodes a "grid()" JSONP callback around the returned json data making it fairly specific to Wax client libraries.

So I went down two paths, first trying to export all the tiles out of mbtiles to json and png files (for those times when you just want to serve static files), then trying to write a simple server that would do dynamic jsonp callbacks. Turns out that in the process, I was able to abstract a lot of the python< ->sqlite interaction into some generic python classes.

Thus python-mbtiles was born. It provides a simple mbtiles web server, a conversion script, and some python classes to work with. No frills, no anything really at this point. More an experiment gone right that might be useful to someone out there in GeoPython land. Enjoy and let me know if you have any ideas!

Average Aspect

Sun 18 March 2012

Ever try to figure out what the average aspect of an area is? i.e.

What direction does this hillside face?

Let's say we want to determine the average elevation of an area based on a raster DEM. Just take the arithmetic mean of all the elevation cells contained in the area - a simple zonal statistics problem.

Turns out that aspect is not quite as straightforward. True, we can easily use gdaldem to create an aspect map.

gdaldem aspect elevation.tif aspect.tif

This gives a raster with values in degrees: 0 is north, 90 is east, 180 is south, etc... but note that 360 is north as well. We're dealing with angular units, not linear units.

For example, take a nearly North facing hillside; the left edge is facing slightly NW (350 degrees) while the right edge faces slighty NE (10 degrees).

The arithmetic mean of the aspect values = (350+350+10+10)/4 = 180°. Due south? That's entirely wrong! It doesn't take into account the angular units. For that we need to create grids representing the sin and cos of the aspect.

Luckily you can use the handy gdal_calc.py utility that comes with recent versions of gdal. This allows you to apply numpy's trigonometric functions to a raster...

gdal_calc.py -A aspect.tif --calc "cos(radians(A))" --format "GTiff" --outfile cos_aspect.tif  
gdal_calc.py -A aspect.tif --calc "sin(radians(A))" --format "GTiff" --outfile sin_aspect.tif

Now we can look at the sum of the cos/sin grid cells for our area and take the arctangent according to this python code

import math
avg_aspect_rad = math.atan2(sum(cos_cells), sum(sin_cells))
avg_aspect_deg = math.degrees(avg_aspect_rad)
print avg_aspect_deg

In our example avg_aspect_deg comes out to an aspect of 0 degrees (due north) which is exactly what we'd expect.

Thanks to Dan Patterson for his forum post which clued me into this approach.

Page 4 / 23