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