Script to generate CONUS raster of various OSM metrics connected to TIGER extent. Counties2012 table is in spherical Albers projection: `+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +nadgrids=@null +units=m +no_defs`

rasterize.py

from itertools import product
from multiprocessing import Pool

from numpy import zeros, ubyte
from psycopg2 import connect
from osgeo import gdal, osr
from PIL import Image

def measure_quad(x, y, cellsize):
    '''
    '''
    global db
    
    if not db:
        db = connect(database='raster_osm', user='osm').cursor()
    
    db.execute('SELECT SRID(geometry) FROM counties2012 LIMIT 1')
    (srid,) = db.fetchone()
    
    ulx, uly = x - 10, y + cellsize + 10
    lrx, lry = x + cellsize + 10, y - 10
    
    bbox = 'SetSRID(MakeBox2D(MakePoint(%(ulx)d, %(uly)d), MakePoint(%(lrx)d, %(lry)d)), %(srid)d)' % locals()
    
    roads = '''('motorway', 'motorway_link', 'trunk', 'trunk_link',
                'primary', 'primary_link', 'secondary', 'secondary_link',
                'tertiary', 'tertiary_link', 'residential', 'road',
                'unclassified')'''
    
    q = '''SELECT Length(Collect(Intersection(Transform(way, %(srid)s), %(bbox)s))), is_touched, osm_uid
           FROM (
             SELECT way, osm_uid,
                    --
                    -- TIGER edits case statement from
                    -- https://github.com/MapQuest/TIGER-Edited-map/blob/master/inc/layer-tiger.xml.inc
                    --
                    (case when osm_uid = '7168' -- DaveHansenTiger
                                and osm_timestamp::timestamp >= '2007-08-03'::timestamp
                                and osm_timestamp::timestamp <= '2008-05-04'::timestamp
                              then 0
                          when osm_uid = '15169' -- Milenko
                                and osm_timestamp::timestamp >= '2007-10-29'::timestamp
                                and osm_timestamp::timestamp <= '2007-12-12'::timestamp
                              then 0
                          when osm_uid = '20587' -- balrog-kun
                                and osm_timestamp::timestamp >= '2010-03-21'::timestamp
                                and osm_timestamp::timestamp <= '2010-04-08'::timestamp
                                and osm_version::int < 3 -- maybe someone else edited between import and name expansion
                              then 0
                          else 1 end) as is_touched
             FROM usa_osm_line
             WHERE way && Transform(%(bbox)s, 4326)
               AND highway IN %(roads)s
           ) AS osm
           GROUP BY is_touched, osm_uid''' \
         % locals()

    db.execute(q)
    
    osm_length, osm_uids = 0.0, 0
    
    for (length, is_touched, osm_uid) in db.fetchall():
    
        osm_length += float(length or 0)
        
        if is_touched and length > 100:
            osm_uids += 1
    
    q = '''SELECT Length(Collect(Intersection(geometry, %(bbox)s)))
           FROM tiger2007
           WHERE geometry && %(bbox)s
             AND ROADFLG = 'Y' ''' \
         % locals()

    db.execute(q)
    row = db.fetchone()
    tiger2007_length = row[0] or 0
    
    q = '''SELECT Length(Collect(Intersection(geometry, %(bbox)s)))
           FROM tiger2012
           WHERE geometry && %(bbox)s
             AND ROADFLG = 'Y' ''' \
         % locals()

    db.execute(q)
    row = db.fetchone()
    tiger2012_length = row[0] or 0
    
    return osm_length, osm_uids, tiger2007_length, tiger2012_length

def make_datasource(filename, array, sref, xform):
    '''
    '''
    driver = gdal.GetDriverByName('GTiff')
    ds = driver.Create(filename, array.shape[1], array.shape[0], 1, gdal.GDT_Float32)
    ds.SetProjection(sref.ExportToWkt())
    ds.SetGeoTransform(xform)
    ds.GetRasterBand(1).WriteArray(array, 0, 0)
    ds.FlushCache()

if __name__ == '__main__':

    db = connect(database='raster_osm', user='osm').cursor()
    
    cellsize = 1000
    
    q = '''SELECT SRID(g), XMin(g), YMin(g), XMax(g), YMax(g)
           FROM (
             SELECT Buffer(Envelope(Collect(geometry)), %(cellsize)d) AS g
             FROM counties2012
          -- WHERE STATEFP = '09'
          -- WHERE STATEFP = '06' AND COUNTYFP = '075'
             WHERE STATEFP NOT IN ('02', '15', '60', '66', '69', '72', '78')
           ) AS env''' \
         % locals()
    
    db.execute(q)
    
    srid, xmin, ymin, xmax, ymax = db.fetchone()
    
    db.execute('SELECT proj4text FROM spatial_ref_sys WHERE srid = %d' % srid)
    
    (proj4,) = db.fetchone()
    
    sref = osr.SpatialReference()
    sref.ImportFromProj4(proj4)
    
    xs = range(int(xmin - (xmin % cellsize)), int(xmax), cellsize)
    ys = range(int(ymin - (ymin % cellsize)), int(ymax), cellsize)
    
    shape = len(ys), len(xs)
    xform = xs[0], cellsize, 0, ys[-1], 0, -cellsize
    
    print 'building arrays...', shape
    
    array_lens = zeros(shape, dtype=float)
    array_uids = zeros(shape, dtype=float)
    array_2007 = zeros(shape, dtype=float)
    array_2012 = zeros(shape, dtype=float)
    
    print 'collecting cells...'
    
    db = None
    
    cells = product(enumerate(xs), enumerate(reversed(ys)))
    
    def multi_quad(((col, x), (row, y))):
        print row, col
        osm_len, osm_uids, tig07_len, tig12_len = measure_quad(x, y, cellsize)
        return row, col, osm_len, osm_uids, tig07_len, tig12_len
    
    results = Pool(8).map(multi_quad, cells, chunksize=10)
    
    print 'populating arrays...'
    
    for (row, col, osm_len, osm_uids, tig07_len, tig12_len) in results:
        array_lens[row, col] += osm_len
        array_uids[row, col] += osm_uids
        array_2007[row, col] += tig07_len
        array_2012[row, col] += tig12_len
    
    print 'saving rasters...'
    
    make_datasource('out-osm-lengths.tif', array_lens, sref, xform)
    make_datasource('out-osm-users.tif', array_uids, sref, xform)
    make_datasource('out-2007-lengths.tif', array_2007, sref, xform)
    make_datasource('out-2012-lengths.tif', array_2012, sref, xform)