As I’m writing this blog about something I did almost 3 years ago, I’m having some trouble remembering the exact circumstances that led me to simulate traffic in a city. However, the general direction was that I was interested in measuring how certain parameters (speed-limit, driver reaction-time, vehicle density etc…) affect congestion (and btw, the definition of congestion is also not trivial).

So, the basic idea is to have some representation of a map - a graph of streets with various metadata, and on this map, vehicles would spawn, and follow some simple rules that model a driver (accelerate to meet the speed limit, decelerate/break when entering a junction or in response to a vehicle in front).

To start, I needed a playing field - a map. So I went to OpenStreetMap a pulled a map file of the Tel-Aviv area.

# Parsing the map

The file is an XML file, and is pretty huge relative to the size of the area it represents. Mainly as it contains a lot of information that’s not relevant to the simulation (street names, buildings etc…).

What I really need (initially) are coordinates of streets. This can be obtained by parsing the <node> tags:

<node id="139698" lat="32.0938800" lon="34.7912317" version="19" timestamp="2011-12-12T17:34:51Z" changeset="10101015" uid="385027" user="Ori952"/>


What we care about are the node’s id, lon (longitude) and lat (latitude). Which can be parsed and stored into nodes like this:

with open('map.osm', 'rU') as xml_file:
tree = etree.parse(xml_file)
root = tree.getroot()
nodes = {}
ways = []
for child in root:
if child.tag == 'node':
lon = float(child.get('lon'))
lat = float(child.get('lat'))
nodes[child.get('id')] = (lon, lat, )


The next type of information we need are the actual streets that connect the nodes. Those are represented by <way> tags:

<way id="5013379" version="16" timestamp="2014-11-15T12:53:09Z" changeset="26797570" uid="331840" user="yrtimiD">
<nd ref="1628415541"/>                                                                                                                                                                       <nd ref="2470788072"/>
<nd ref="2109540860"/>
<tag k="highway" v="pedestrian"/>
<tag k="int_name" v="Nahalat Binyamin"/>
<tag k="motor_vehicle" v="no"/>
<tag k="name" v="נחלת בנימין"/>
<tag k="name:en" v="Nachalat Biniamin"/>
<tag k="name:he" v="נחלת בנימין"/>
</way>


Initially, what we need is just the list of node ids that compose the road. However, there are two important points:

• Not all <way>s represent vehicle roads. This is determined by the v (value) property of the <tag> whose k (key) is highway:

NON_ROAD_TYPES = ['proposed', 'platform', 'path', 'footway', 'pedestrian',
'cycleway', 'road', 'steps', 'track', 'unclassified', 'construction']
for child in way:
if child.tag == 'tag' and child.get('k') =='highway':
return False
return True
return False

• Some roads are two-way. As our graph will be directed, we need to manually add the arcs in the opposite direction.

Eventually we can parse the road data into the ways list using:

for child in root:
if child.tag == 'way':
way = child
continue
nds = []
for nd in way:
if nd.tag == 'nd':
nds.append(nd.get('ref'))
ways.append((nds, speed))
ways.append((list(reversed(nds)), speed))


Now, we’re almost there, as we have a list of “ways” which are themselves lists of nodes. But what we actually need is arcs, which are tuples of consecutive nodes:

ref_arcs = []
for way in ways:
arc_tuples = zip(way[:-1], way[1:])
ref_arcs.extend((src, dst) for src, dst in arc_tuples)


# Refining

If you think about it, the map file we use is a cropped map. So what would happen is that we will have lots of roads that leave or exit the map that are dead-end. If we don’t want our map to have car sinks, we need to eliminate all those dead-ends:

while True:
src_set = set(src for src, _, _ in ref_arcs)
_ref_arcs = []
for src, dst, speed in ref_arcs:
if dst in src_set:
_ref_arcs.append((src, dst, speed))
else:
# If no arcs were discarded in this round we are done
break
ref_arcs = _ref_arcs


# Coordinates

The node coordinates are all angular - longitude and latitude. However, on the ground, vehicles are moving at km/h (or m/s). We have no way of running the simulation without assigning fixed coordinated to the nodes.

There is an entire science to measuring the earth. But I wanted to choose something simple.

Basically what I did was take the lon/lat bounding box the of map. This describes a part of the earth’s shell:

Now, we draw a plane tangent to the earth at the center of that “shell”:

Now, suppose we have a point $\vec E$ on the earth’s surface:

We will cast a ray from the center of the earth to that point, and continue tracing it until it intersects with the plane (or map if you will) at $\vec M$ :

We now have a 1-to-1 mapping from points on the surface of the earth to points on a flat map!

Now we need a way to get 2D coordinates. To do that we need to select axes.

We start with a north “hint”, we can use the $\hat z$ vector as it points at the north pole. Obviously, this vector does not fit in map, but sticks out of it. We want to take the component of $\hat z$ that’s embedded in the map and use that as “north”. Generally, a vector can be decomposed in relation to a plane to a perpendicular (in the direction of the plane’s normal) and transverse (embedded in the plane, and thus perpendicular to the normal) $\hat z$ :

Since $\vec z_{\parallel}$ is perpendicular to $\hat n$ , then $\vec z_{\parallel} \cdot \hat n = 0$ .

In addition, $\vec z_{\bot}$ is in the direction of $\hat n$ , so we can write it as: $\vec z_{\bot} = \alpha \hat n$

We can now use all that to extract $\vec z_{\parallel}$ (which will serve as the plane’s north):

On the surface of the map, we can have east/north axis and thus a 2D coordinate system.

So, first, I want to convert lon/lat coordinates to a 3d point on the surface of earth’s sphere. These are just Spherical Coordinates:

def vector3_from_lon_lat(lon, lat):
l = lon * math.pi / 180
p = lat * math.pi / 180

return (math.cos(p) * math.cos(l), math.cos(p) * math.sin(l), math.sin(p))


What this would give us is the normal to the tangent plane: $\hat n_0 = \hat r(\frac{\pi}{180}lon, \frac{\pi}{180}lat)$. Which will give us the following equation for the tangent plane:

# Bounding box
minlon = min(lon for lon, lat in nodes.values())
maxlon = max(lon for lon, lat in nodes.values())
minlat = min(lat for lon, lat in nodes.values())
maxlat = max(lat for lon, lat in nodes.values())

# Center point
lon_0 = (minlon + maxlon) / 2
lat_0 = (minlat + maxlat) / 2

n_0 = vector3_from_lon_lat(lon_0, lat_0)


Now we can take the $\hat z$ and project it onto the tangent plane to get a “north” vector:

z = (0, 0, 1) # "real" north: points to the north pole

north_0 = sub3(z, scale3(n_0, dot3(z, n_0))) # Project north onto the tangent plane
north_0 = scale3(north_0, 1 / norm3(north_0)) # Normalize


East is then obtained by using a cross product (right-hand rule):

east_0 = cross3(north_0, n_0) # Right-hand rule


The north and east facing unit-vectors on the plane will be used as basis vectors.

Now for each lon/lat coordinate, we scale it using earth’s radius (6.317e6[m]) towards the plane:

If we plug this into the equation for the tangent plane:

Which gives us the following intersection point on the plane:

We can then decompose to north and south coordinates (in meters!):

n = vector3_from_lon_lat(lon, lat)
t = R_E / dot3(n, n_0)
x = scale3(n, t)
u = dot3(x, north_0)
v = dot3(x, east_0)


That’s it so far, next time we’ll talk about setting up a simulation.