Elijah Robison | GIS Blog

A scrapbook of GIS tricks, with emphasis on FOSS4G.

Archive for the ‘Map Tiling’ Category

Create UTFGrid Tiles from PostGIS Tables

without comments

I assume you’re where I was about a week ago. That is, you’ve heard of UTFGrid, and now you want to render your own UTFGrid tiles. Perhaps you found yourself looking at this thread over at GIS.SE, but you don’t want to jump into TileStache just now. Anyway if you’ve got a working installation of GDAL/OGR and Mapnik 2+, complete with Python bindings, I’ll show you what worked for me..

Because this is merely an adaptation of Matthew Perry’s original solution, I highly recommend considering his original blog entry on the topic, complete with discussion points and caveats, before proceeding further!

Once you’re ready, head over to GitHub and download some code, specifically globalmaptiles.py, create_tile_shp.py, and create_utfgrids.py. [1] You don’t actually need create_tile_shp.py, but I snagged all three. To keep things simple put these files in the same folder.

Next, in that same directory, create a new Python file, I called mine createUtfgridsFromPG.py.

Where create_utfgrids.py works entirely on shapefiles, createUtfgridsFromPG.py accepts an OGR PostgreSQL connection string in place of a shapefile path. Fortunately the original OGR code didn’t change, but to use the Mapnik PostGIS driver I had to iterate over the PostgreSQL connection string and store the connection parameters so I could supply them a differently to Mapnik.

Finally, copy the following code and paste it into your new file. It’s basically the same as create_utfgrids.py, but I changed shppath to pgconn and added some code to use a mapnik.PostGIS() datasource in place of a mapnik.Shapefile() datasource.

If you’re curious, the comment # ELR 2014.9.26: flags the few places I changed the original code.

createUtfgridsFromPG.py

#!/usr/bin/env python
# -*- coding: utf-8  -*-
"""
create_utfgrids.py
Author: Matthew Perry
License: BSD

Creates utfgrid .json tiles for the given polygon shapefile

Thx to Dane Springmeyer for the utfgrid spec and mapnik rendering code
and to  Klokan Petr Přidal for his MapTiler code
(http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection/)

"""
import globalmaptiles
import mapnik
import ogr
import os
from optparse import OptionParser, OptionError
try:
    import simplejson as json
except ImportError:
    import json

def create_utfgrids(pgconn, minzoom, maxzoom, outdir, fields=None, layernum=0):

    # ELR 2014.9.26:
    # Original implementation pushed in a shapefile path.
    #ds = ogr.Open(shppath)
    ds = ogr.Open(pgconn)

    # ELR 2014.9.26:
    # Iterate over the PostgreSQL connection string and pull out values we need
    # to use Mapnik's PostGIS datasource constructor.
    pgConnARR = pgconn[3:].split(' ')
    for kvPair in pgConnARR:
        if kvPair.split('=')[0] == "host":
            nikHost = kvPair.split('=')[1]
        if kvPair.split('=')[0] == "port":
            nikPort = kvPair.split('=')[1]
        if kvPair.split('=')[0] == "user":
            nikUser = kvPair.split('=')[1]
        if kvPair.split('=')[0] == "password":
            nikPass = kvPair.split('=')[1]
        if kvPair.split('=')[0] == "dbname":
            nikDB = kvPair.split('=')[1]
        if kvPair.split('=')[0] == "tables":
            nikTable = kvPair.split('=')[1]

    print
    print "WARNING:"
    print " This script assumes a polygon shapefile in spherical mercator projection."
    print " If any of these assumptions are not true, don't count on the results!"
    # TODO confirm polygons
    # TODO confirm mercator
    # TODO get layernum from command line
    layer = ds.GetLayer(layernum)
    bbox = layer.GetExtent()
    print ""
    print str(bbox)

    mercator = globalmaptiles.GlobalMercator()

    m = mapnik.Map(256,256)

    # Since grids are `rendered` they need a style
    s = mapnik.Style()
    r = mapnik.Rule()
    polygon_symbolizer = mapnik.PolygonSymbolizer(mapnik.Color('#f2eff9'))
    r.symbols.append(polygon_symbolizer)
    line_symbolizer = mapnik.LineSymbolizer(mapnik.Color('rgb(50%,50%,50%)'),0.1)
    r.symbols.append(line_symbolizer)
    s.rules.append(r)
    m.append_style('My Style',s)

    print ""
    # ELR 2014.9.26:
    # Original implementation using shapefile..
    #ds = mapnik.Shapefile(file=shppath)

    # ELR 2014.9.26:
    # Parameterized PostGIS implementation..
    ds = mapnik.PostGIS(host=nikHost,port=nikPort,user=nikUser,password=nikPass,dbname=nikDB,table=nikTable)

    mlayer = mapnik.Layer('poly')
    mlayer.datasource = ds
    mlayer.styles.append('My Style')
    m.layers.append(mlayer)

    print ""
    if fields is None:
        fields = mlayer.datasource.fields()
        print "Fields were NONE. Using.."
        print fields
    else:
        print "Fields are USER PROVIDED. Using.."
        print fields
    print ""

    for tz in range(minzoom, maxzoom+1):
        print " * Processing Zoom Level %s" % tz
        tminx, tminy = mercator.MetersToTile( bbox[0], bbox[2], tz)
        tmaxx, tmaxy = mercator.MetersToTile( bbox[1], bbox[3], tz)
        for ty in range(tminy, tmaxy+1):
            for tx in range(tminx, tmaxx+1):
                output = os.path.join(outdir, str(tz), str(tx))
                if not os.path.exists(output):
                    os.makedirs(output)

                # Use top origin tile scheme (like OSM or GMaps)
                # TODO support option for TMS bottom origin scheme (ie opt to not invert)
                ymax = 1 << tz;
                invert_ty = ymax - ty - 1;

                tilefilename = os.path.join(output, "%s.json" % invert_ty) # ty for TMS bottom origin
                tilebounds = mercator.TileBounds( tx, ty, tz)
                #print tilefilename, tilebounds

                box = mapnik.Box2d(*tilebounds)
                m.zoom_to_box(box)
                grid = mapnik.Grid(m.width,m.height)
                mapnik.render_layer(m,grid,layer=0,fields=fields)
                utfgrid = grid.encode('utf',resolution=4)
                with open(tilefilename, 'w') as file:
                    file.write(json.dumps(utfgrid))

if __name__ == "__main__":
    usage = "usage: %prog [options] shapefile minzoom maxzoom output_directory"
    parser = OptionParser(usage)
    parser.add_option("-f", '--fields', dest="fields", help="Comma-seperated list of fields; default is all")
    (options, args) = parser.parse_args()

    if len(args) != 4:
        parser.error("Incorrect number of arguments")

    pgconn = args[0]
    minzoom, maxzoom = int(args[1]), int(args[2])
    outdir = args[3]

    if os.path.exists(outdir):
        parser.error("output directory exists already")

    if options.fields:
        fields = options.fields.split(",")
    else:
        fields = None

    create_utfgrids(pgconn, minzoom, maxzoom, outdir, fields)

Usage..

Once you’ve prepared createUtfgridsFromPG.py, you can call it from the command line like this..

C:\xDev\utfgrids\createUtfgridsFromPG.py "PG:host=127.0.0.1 port=5432 user=postgres dbname=gis password=passw0rd tables=parcels_pmerc" 12 16 "C:/xGIS/tiles/utf" -f tms,owner_name

  • Hopefully the PostgreSQL connection string ("PG:host=..") makes sense.
  • 12 and 16 represent the minimum and maximum zoom levels to be rendered, respectively.
  • The directory “C:/xGIS/tiles/utf” is where your UTFGrid tiles will be saved.
  • And -f tms,owner_name,the_wkt represents a comma-separated list of data fields you want in your UTFGrid.

Caveats..

  • Both create_utfgrids.py and createUtfgridsFromPG.py require your geodata table to be in a Web Mercator projection (EPSG:3857)!
  • The script assumes a top-origin tile scheme, like OSM and others.
  • The script will only work with polygons.
  • While the OGR PostgreSQL connection string has a tables parameter, this implementation will only accept one table.
  • The script will create your target directory, in the example case, utf, and it will throw an error if you create this directory in advance.

[1] Many thanks to Matthew Perry, Klokan Petr Přidal, and Dane Springmeyer for their collective efforts and for sharing their work.

Written by elrobis

September 26th, 2014 at 1:26 pm