Geospatial Representation of Roads

It turns out that the polygons in the Roads shapefiles do not represent roads the way I thought they would.  To explore the shapefile, I wrote a script using pyshp: explore_shapefile.py.  The pyshp module is definitely more rewarding to use than GDAL, since it is pure Python and it is thus possible to use Python's type introspection to explore.  I've found that all of the scripts I've written in Python 3 will run also in Python 2, which is convenient.  I also came across one of the big improvements Python 3 has over Python 2 in terms of type introspection.

The following line of code prints out the type of the first element of a list from a shapefile: 

print('type(sr[0]): {0}'.format(type(sr[0])))

Running the script in Python 2 results in: 

type(sr[0].shape): <type 'instance'>

which is not very informative.  Python 3 does much better: 

type(sr[0]): <class 'shapefile._ShapeRecord'>

Anyway, since my final project is due tomorrow, I had to cut the investigation into using pyshp short, but I'll be sure to continue looking into it soon.

My goal was to represent Columbia Pike in bright yellow or orange on my map. Columbia Pike is the road on which the proposed streetcar that had such a profound impact on the Vihstadt / Howze election ran.  I thought adding it to the map would provide visual interest.  It also gave me a chance to explore adding more than one layer to a map.

I wrote several scripts to make maps as part of investigating this.

draw_roads.py

Renders all 883 shapes from the Roads shapefile in a map.

a_road_shape.py

Renders a single shape from the Roads shapefile as a new layer on the precincts map.

roads_near_columbia_pike.py

I used trial and error (and a bit of divide and conquer) to find shapes among the 883 that represented roads on and near Columbia Pike.

election_map_with_roads_near_cp.py

Finally, I added this roads layer to the election results map I made yesterday.

The first script, draw_roads.py, generates a map of the roads in Arlington County:

I have yet to begin studying QGIS, but I googled how to load a shapefile as a vector layer to see how it would compare:

Since I generally feel more comfortable expressing my thoughts in code rather than using tools, I think I'm going to generally prefer writing python scripts to using qgis to make maps, but I realize I'll have to do both.

I assumed that each of the 883 shapes in the Roads shapefile represented a single road.  I was disappointed that the metadata didn't include names associated with each of the shapes, but I figured I could use a kind of binary search to print out half of them, render the map, and see which half Columbia Pike was in, cut that in half, and so forth, until I found Columbia Pike.  Soon I'll be able to read geospatial coordinates from the shapes and write a script that selects the shapes I want, but for now I figured this approach would be the quickest way to get to my goal.

It turns out that the reason there are no road names associated with the shapes is that each shape represents a cluster of connected road segments rather then a single road.  a_road_shape.py illustrates this:

The orange line segments are from a single shape from the Roads shapefile.  Once I realized that, I spent some time gathering shapes that had a piece of the Pike on them, and made roads_near_columbia_pike.py to render them:

Finally, I merged this with the color coded election results map I made yesterday in election_map_with_roads_near_cp.py to produce:

The color contrast is bad here, but I'm running out of time to complete the project, and since this fails to render just Columbia Pike the way I wanted anyway, I decided not to spend any more time on it.

Where I'll Go From Here

OpenStreetMap is an incredible resource for open source mapping projects, and yesterday I did a quick search to see if it might provide me with a solution to my rendering Columbia Pike project.  I found out that OSM stores streets as a collections of ordered nodes called ways. This way holds a piece of Columbia Pike near where I live:

https://www.openstreetmap.org/way/334620323#map=18/38.85606/-77.10834

Here are the things I think I should investigate next:

• OSM data
• pyproj

My immediate project goal will be to render Columbia Pike, this time from OSM, onto the election map as I had originally intended.  I anticipate that the pyproj module will be the tool to use to convert among different geospatial data formats.