Red Street, Blue Street

In the end, I went with the upgrade to Mathematica 7. Of all the new features, the one that really hooked me–which is comparatively minor, compared to all the other new features–is the ability to import SHP files. The importation is not terribly well documented nor is there much additional support, but it was pretty easy to do a few nifty things with the DC Street Centerline file.

As you may know, there is a street in DC for every state in the union. Pennsylvania Avenue is probably the most famous of these; the White House sits at 1600 Pennsylvania Ave NW. I used to live on Massachusetts Avenue. So my first idea was to make a street map of DC in which the state-named streets were colored red-ish or blue-ish depending on their vote in the recent election.

Here it is:

Read on to see how I made it:

I imported election returns from Wikipedia, but the process of picking the data I wanted out of the whole Import was so ugly that I’m not going to step through it here. I was inspired, though, by this similar post on the Wolfram blog. If anyone at Wolfram is listening: it’d be really great if we could import data from the web by, say, clicking on a table of a web page displayed in a browser, instead of picking out which item in some nested list is the one we want.

From the raw vote data, I computed the (approximate) percentage of the vote that Obama got:¹

ntrs[x_] := N[ToExpression[StringReplace[x, "," -> ""]]]
percentData = {ToUpperCase[#[[1]]], ntrs[#[[2]]]/(ntrs[#[[2]]] + ntrs[#[[3]]])} & /@ votedata

Which gives a list of things like

{"CALIFORNIA", 0.622799}.
The conversion to uppercase is to match the way the DC street names are recorded.

I’ll need a list of states, in order to match the street names. I get this from Mathematica’s curated CountryData data set, using a Regular Expression to insert a space any time there’s a lowercase letter followed by an uppercase letter, to change, for example, “NewJersey” into “New Jersey”:

states = StringReplace[#, RegularExpression["([a-z])([A-Z])"] -> "$1 $2"] & /@ CountryData["United States", "Regions"]


And the street datafile itself. When Mathematica imports SHP files, it returns a list of layers; with only one layer in this file, we get a (long) list with one element. So I specify in the Import that I only want this first (and only) element of the list:

dcstreets=Import["StreetSegmentsLn.shp","Data"][[1]];

And this gives a list of 4 replacement rules, some of which are very long:

In[]:=dcstreets[[1]]
Out[]:=LayerName -> StreetSegmentLn
In[]:=dcstreets[[2]]
Out[]:=Geometry -> long list of graphics primitives defining the data elements
In[]:=dcstreets[[3]]
Out[]:=Labels -> list of the names of metadata fields
In[]:=dcstreets[[4]]
Out[]:=LabeledData -> list of (long) lists of metadata for each data element

For this data file, the data elements are segments of streets between intersections, so the graphics primitives are lines defined by two (or more, for curved streets) points. There are 25 metadata fields, including the name of each street, the address range covered by the segment, the quadrant, and the type of road (e.g. Local, Principal Arterial, etc.). The FullForm of Mathematica replacement rule, such as lhs -> rhs, is Rule[lhs,rhs], so we access the data in each of the four rules as the second element of another list. Thus we can get at, say, the graphics primitives for the 27th street segment as

In[]:=dcstreets[[2,2,27]]
Out[]:=Line[{{396233., 138045.}, {396233., 138132.}, {396233., 138142.}}]

For ease in working with these data, I combined the geometry and metadata into one table for the 13414 street segments:

dcstreetdata = Transpose[Join[{dcstreets[[2, 2]]}, dcstreets[[4, 2, All, 2]]]];

I separate this into two groups: state-named streets, and the rest:

statestreets = Select[dcstreetdata, MemberQ[ToUpperCase[states], #[[6]]] &];
nonstatestreets = Select[dcstreetdata, ! MemberQ[ToUpperCase[states], #[[6]]] &];

The 6th element of each entry in dcstreetdata is the name of the street. From this, I create two graphics:

nonstatemap = Graphics[{Gray, nonstatestreets[[All, 1]]}]


statemap = Graphics[Function[x, Tooltip[{Directive[Thick, RGBColor[.5 + 2 (.5 - x[[2]]), 0, .5 + 2 (x[[2]] - .5)]], Transpose[Select[statestreets, Function[y, y[[6]] == x[[1]]]]][[1]]}, x[[1]]]] /@ percentData]

I’ve packed a lot into statemap, using nested pure functions: it makes a Tooltip with the state/street name popping up upon mouseover (it works better in Mathematica than in the export HTML). For each state from the list of election returns in percentData, it finds all the street segments named after that state and makes a graphic element out of the geometry of that street segment, of a color that ranges from pure blue to pure red. The red-blue divide is stretched so that pure blue corresponds to the percentage of Obama’s highest percentage (Hawaii, 72.99%) and pure red corresponds to Obama’s lowest percentage (Wyoming, 33.43%). The final map is created with Show:

votemap = Show[statemap, nonstatemap]

So far, so good. I’ll write more on Mathematica and DC GIS data soon.

1. really, it’s Obama’s percentage of the combined Obama/ McCain vote, ignoring votes for other candidates [↩]