Showing posts with label visualization. Show all posts
Showing posts with label visualization. Show all posts

Tuesday, October 5, 2010

Beta release: How to render OpenStreetMaps using Python and Matplotlib

Introduction:

Open Street Maps is an collaboratively developed and open source map database. The entire OSM database or specific portions can be downloaded. The database includes information on roads, political boundaries, natural features, rail lines, trails, and much more. The website can generate output files with sufficient information to create custom and detailed surface street maps for many areas of the world. With Python and Matplotlib the XML data in an OpenStreetMaps OSM file can be rendered as a custom map can be rendered to meet specific needs. Below is an example of a map which can be generated from and OSM data file.


This is a small snippet of code which makes it relatively easy to see how to use the raw data from Open Street Maps. For a complete rendering examples see the Open Street Maps wiki and the pyrender source code.

Key Concepts:

An OSM file is generated by selecting an area of interest at www.OpenStreetMap.org and selecting the export tab. Under the export tab, there is an option to export the map to 'OpenStreetMap XML Data'. If this is selected, then an XML file which contains all of the data for the selected region is generated. This XML file includes road, railways, trails, and land features among other things. The XML file has a a few header elements which define the XML standard, the data source, and the bounds on the data set. After that, the data file consists of a long list of nodes in the map. Each node describes one point on the surface of the earth. After the nodes, the ways are defined. A way consists of a list of nodes and tags. The nodes are are order and describe something like a road or a lake. The tag provide sufficient data to understand what the way represents. After the ways, come relations. Relations combine other elements and provide tags which area associated with those elements.
To render roads, only nodes and ways need to be considered.
When plotting lots of elements, Matplotlib can be sped up significantly (3X or more) by turning off autoscaling.
Since the ways are not ordered in the OSM file, z-ordering in matplotlib provides a nice way to control what has the most prominence.
So that the roads which are wider than 1 pixel have smooth transitions, the solid_capstyle and the solid_endstyle are set to ‘round’.
To simplify conversion of the XML data in the OSM file into an object which behaves like a dictionary, a script from ActiveState was used.

The Code:

download source code
download example OSM file
#############################################################
## from http://code.activestate.com/recipes/534109-xml-to-python-data-structure/

import re
import xml.sax.handler

def xml2obj(src):
    """
    A simple function to converts XML data into native Python object.
    """

    non_id_char = re.compile('[^_0-9a-zA-Z]')
    def _name_mangle(name):
        return non_id_char.sub('_', name)

    class DataNode(object):
        def __init__(self):
            self._attrs = {}    # XML attributes and child elements
            self.data = None    # child text data
        def __len__(self):
            # treat single element as a list of 1
            return 1
        def __getitem__(self, key):
            if isinstance(key, basestring):
                return self._attrs.get(key,None)
            else:
                return [self][key]
        def __contains__(self, name):
            return self._attrs.has_key(name)
        def __nonzero__(self):
            return bool(self._attrs or self.data)
        def __getattr__(self, name):
            if name.startswith('__'):
                # need to do this for Python special methods???
                raise AttributeError(name)
            return self._attrs.get(name,None)
        def _add_xml_attr(self, name, value):
            if name in self._attrs:
                # multiple attribute of the same name are represented by a list
                children = self._attrs[name]
                if not isinstance(children, list):
                    children = [children]
                    self._attrs[name] = children
                children.append(value)
            else:
                self._attrs[name] = value
        def __str__(self):
            return self.data or ''
        def __repr__(self):
            items = sorted(self._attrs.items())
            if self.data:
                items.append(('data', self.data))
            return u'{%s}' % ', '.join([u'%s:%s' % (k,repr(v)) for k,v in items])

    class TreeBuilder(xml.sax.handler.ContentHandler):
        def __init__(self):
            self.stack = []
            self.root = DataNode()
            self.current = self.root
            self.text_parts = []
        def startElement(self, name, attrs):
            self.stack.append((self.current, self.text_parts))
            self.current = DataNode()
            self.text_parts = []
            # xml attributes --> python attributes
            for k, v in attrs.items():
                self.current._add_xml_attr(_name_mangle(k), v)
        def endElement(self, name):
            text = ''.join(self.text_parts).strip()
            if text:
                self.current.data = text
            if self.current._attrs:
                obj = self.current
            else:
                # a text only node is simply represented by the string
                obj = text or ''
            self.current, self.text_parts = self.stack.pop()
            self.current._add_xml_attr(_name_mangle(name), obj)
        def characters(self, content):
            self.text_parts.append(content)

    builder = TreeBuilder()
    if isinstance(src,basestring):
        xml.sax.parseString(src, builder)
    else:
        xml.sax.parse(src, builder)
    return builder.root._attrs.values()[0]



#############################################################



def render(myMap):



    # make dictionary of node IDs
    nodes = {}
    for node in myMap['node']:
        nodes[node['id']] = node
        
    ways = {}
    for way in myMap['way']:
        ways[way['id']]=way



    import pylab as p


    renderingRules = {
        'primary': dict(
                linestyle       = '-',
                linewidth       = 6,
                color           ='#ee82ee', 
                zorder          = -1,
                ),
        'primary_link': dict(
                linestyle       = '-',
                linewidth       = 6,
                color           = '#da70d6',
                zorder          = -1,            
                ),
        'secondary': dict(
                linestyle       = '-',
                linewidth       = 6,
                color           = '#d8bfd8',
                zorder          = -2,            
                ),
        'secondary_link': dict(
                linestyle       = '-',
                linewidth       = 6,
                color           = '#d8bfd8',
                zorder          = -2,            
                ),
        'tertiary': dict(
                linestyle       = '-',
                linewidth       = 4,
                color           = (0.0,0.0,0.7),
                zorder          = -3,            
                ),
        'tertiary_link': dict(
                linestyle       = '-',
                linewidth       = 4,
                color           = (0.0,0.0,0.7),
                zorder          = -3,            
                ),
        'residential': dict(
                linestyle       = '-',
                linewidth       = 1,
                color           = (0.1,0.1,0.1),
                zorder          = -99,            
                ),            
        'unclassified': dict(
                linestyle       = ':',
                linewidth       = 1,
                color           = (0.5,0.5,0.5),
                zorder          = -1,            
                ),
        'default': dict(
                linestyle       = '-',
                linewidth       = 3,
                color           = 'b',
                zorder          = -1,            
                ),
        }
                        

    # get bounds from OSM data            
    minX = float(myMap['bounds']['minlon'])
    maxX = float(myMap['bounds']['maxlon'])
    minY = float(myMap['bounds']['minlat'])
    maxY = float(myMap['bounds']['maxlat'])



    fig = p.figure()
    
    # by setting limits before hand, plotting is about 3 times faster
    ax = fig.add_subplot(111,autoscale_on=False,xlim=(minX,maxX),ylim=(minY,maxY))
    
    for idx,nodeID in enumerate(ways.keys()):
        wayTags         = ways[nodeID]['tag']
        if not wayTags==None:
            hwyTypeList  = [d['v'] for d in wayTags if d['k']=='highway']
            if len(hwyTypeList)>0:
                    wayType = hwyTypeList[0]  
            else:
                    wayType = None
        else:
            wayType = None
        try:
            if wayType in ['primary','primary_link',
                            'unclassified',
                            'secondary','secondary_link',
                            'tertiary','tertiary_link',
                            'residential',
                            'trunk','trunk_link',
                            'motorway','motorway_link',
                            ]:
                oldX = None
                oldY = None
                
                if wayType in renderingRules.keys():
                    thisRendering = renderingRules[wayType]
                else:
                    thisRendering = renderingRules['default']
                    
                for nCnt,nID in enumerate(ways[nodeID]['nd']):
                    y = float(nodes[nID['ref']]['lat'])
                    x = float(nodes[nID['ref']]['lon'])
                    if oldX == None:
                        pass
                    else:
                        p.plot([oldX,x],[oldY,y],
                            marker          = '',
                            linestyle       = thisRendering['linestyle'],
                            linewidth       = thisRendering['linewidth'],
                            color           = thisRendering['color'],
                            solid_capstyle  = 'round',
                            solid_joinstyle = 'round',
                            zorder          = thisRendering['zorder'],
                            )
                    oldX = x
                    oldY = y
                       
        except KeyError:
            pass
          
    print'Done Plotting'
    p.show()




########################################

src = file('map.osm')
myMap = xml2obj(src)
render(myMap)

Testing Configuration:

Other Useful Links and references:


Questions Answered:

  • How to visualize an Open Street Map from data
  • How to plot a map
  • How to draw an Open Street Map from the OSM file

Monday, April 6, 2009

Build the working directory

To replicate the work I'm going to share, you will need to build a directory that will share the source code and the data. Follow these steps to build that directory.
  1. Create a new directory for testing.
  2. Download the 2001 NHTS ASCII.CSV data set.
  3. Unzip the contents of that file into the testing directory. You should have four CSV files.

If you do not have Python already available on your machine, I'd recommend starting with Portable Python, I'm using Portable Python 1.1 which implement Python 2.5.4 as a portable program (try saying that 3 time fast!). Everything initially will work in this implementation.

Using Python to Visualize the 2001 NHTS - Project Scope

This project will work through using Python for visualizing a complex data set to gain insight and understanding. One of my favorite data sets is the US National Household Travel Survey. This survey summarizes demographics and a days worth of driving for 70,000 people over a one year period. The data set is so complex that after about two years of looking at it, I've only just started to find the interesting information contained in it. My initial goals are as follows:
  1. Write pure Python libraries for inputting the data sets and preprocessing the results. I'm going to write everything in pure Python from scratch, rather than use existing tools for two reasons. First I need the exercises to learn Python better. Second, I want my tools for importing the data to reusable under Jython and for use with Blender out of the box.
  2. Use the following tools to visualize information in the NHTS:
  3. I'll be using the following Python distributions:
  4. Create both static and animated visualizations in each tool.
  5. Gain some insight into travel patterns that I did not have before starting this project.
  6. Build some tools for exploring the 2008 release of the NHTS.