I have a shapefile of ~8000 polyline vector features representing traffic journeys on road routes. Each has an attribute representing volume of traffic for that journey. The routes are all generated by the same routing engine so the node coordinates should match exactly where routes follow common sections of road. I am wanting to effectively add these together in a "map algebra" sort of way to give a representation of network usage. Keeping the result in vector would be preferable, but ultimately not essential.

I can see how this could be done conceptually simply by converting each feature to a separate raster of pixel values = traffic volume, then summing all the rasters but this would practically be very painful to combine 8000 high-resolution rasters. Ironically, QGIS is effectively already doing this visually for me - I am displaying the layer using combining of features by "addition" method.. but I can't see any way of getting that directly out into a much higher resolution raster format than my display:

QGIS-rendered display screengrab

Surely there must be a vector way to do this? (QGIS or ArcGIS) If it makes it easier to not go the shapefile route I have all the route nodes in CSV format, one line per route. The polygon "union" then "dissolve" method outlined in How to union polygons and add attribute values of combined output features looked promising, but Union does not work on line features it seems. Flicking through the ArcGIS Network Analyst "Route Analysis" pages looked potentially promising if rather OTT for this purpose and I'm guessing would require a lot of messing round to get there - I'm sure there must be a simpler way!


Here is an example of how to perform this using Python, Fiona, and Shapely (I tried doing it in arcpy but gave up in frustration). I'm not sure if this is the most efficient approach but it should do the trick. The function first builds a set of unique segments from all of the routes, then identifies which route each segment is contained by (summing up the volume on each route), and finally writes out the segment and total volume to a new shapefile. The script is a bit rough and only tested on a sample data set I constructed. You will need to modify it to work with your particular data, and I'm not sure how quickly it will run with ~8000 input routes.

If you end up using please post a comment and let me know how long it took to run.

import fiona
from shapely.geometry import shape, LineString, mapping

def explode():
    merged_lines = None
    unique_segments = []
    input_routes = []

    #Read input file and 
    #build set of unique segments across all routes

    with fiona.open('testarcs.shp', 'r') as source:
        source_crs = source.crs
        source_driver = source.driver

    for route in source:
            #cache the current route and get current route geometry
            route_geom = shape(route['geometry'])

            #break current route into segments
            route_segments = [LineString([route_geom.coords[seg_indx], route_geom.coords[seg_indx + 1]]) for seg_indx in xrange(len(route_geom.coords) - 1)]

            if not merged_lines:
                #this is the first iteration so all segments are unique
                merged_lines = route_geom
                #test segments to see if they are coincident with the current merged routes
                unique_segments.extend([curr_seg for curr_seg in route_segments if not merged_lines.contains(curr_seg)])
                #union the current route with all routes that have been processed so far
                merged_lines = merged_lines.union(route_geom)

    output_schema = {'geometry':'LineString','properties':{'Tot_Vol':'int'}}

    #create output file
    with fiona.open('seg_summ.shp','w',driver=source_driver,crs=source_crs, schema=output_schema) as seg_outputs:
        new_rec = {}
        #Now iterate through the unique segments.
        #For each segment find the routes that it intersects with and
        #sum the volume on the routes
        for curr_segment in unique_segments:
            curr_volume = 0
            for curr_route in input_routes:
                if curr_route[0].contains(curr_segment):
                    curr_volume += curr_route[1]['properties']['Volume']    

            new_rec['properties'] = {'Tot_Vol' : curr_volume}
            new_rec['geometry'] = mapping(curr_segment) 

def explode_rtree():
  segments = {}
  segments_index = rtree.index.Index()
  segment_id = 0
  source_crs = None
  source_driver = None

  #read in routes and break down to segments
  with fiona.open('testarcs.shp', 'r') as source:
    source_crs = source.crs
    source_driver = source.driver

    for route in source:
      #get current route geometry and volume
      route_geom = shape(route['geometry'])
      route_volume = route['properties']['Volume']

      route_segments = [LineString([route_geom.coords[seg_indx],route_geom.coords[seg_indx + 1]]) for seg_indx in xrange(len(route_geom.coords) - 1)]
      for curr_segment in route_segments:
        #get set of segments with overlapping bounds
        candidates = list(segments_index.intersection(curr_segment.bounds))
        if not candidates:
          #current segment does not intersect with any existing segment bounds
          #this means its a new segment
          segments[segment_id] = [curr_segment,route_volume]
          segment_id += 1
          #check each candidate to see if it the same as the current
          for curr_cand in candidates:
            if segments[curr_cand][0].equals(curr_segment):
              #the current segment is coincident with a segment that
              #has already been processed so add volume
              segments[curr_cand][1] += route_volume
          else: #else clause of for loop
            #the current segment does not intersect with any candidates
            #so add it as a new segment
            segments[segment_id] = [curr_segment,route_volume]
            segment_id += 1

  output_schema = {'geometry':'LineString','properties':{'Tot_Vol':'int'}}

  #create output file
  with fiona.open('seg_summ2.shp','w',driver=source_driver,crs=source_crs,schema=output_schema) as seg_outputs:
    for segment_id in segments:
      new_rec = {}

      new_rec['properties'] = {'Tot_Vol' : segments[segment_id][1]}
      new_rec['geometry'] = mapping(segments[segment_id][0]) 

  • Thanks - got it working (in Windows XP) with result exactly as desired for 10 routes. I had to indent the code for the loop "for route in source" as otherwise the previous "with" line caused the collection source to be closed on completion. Unfortunately the merge method seems too slow - 15 minutes to read in first 50 routes, 30 mins for next 50, 52 mins for next 50. Not going to manage 8000! I imagine generating a hash for each segment could speed up searches significantly, or maybe the nature of the search is inherently slow in a vector format and I need a raster solution? – Bristle6 Aug 4 '14 at 13:16
  • Another option would be to use an rtree index, which would make searching for intersecting routes more efficient. I have not done this yet with shapely and fiona but here's a link that describes how it could be done. snorf.net/blog/2014/05/12/using-rtree-spatial-indexing-with-ogr – dblanchett Aug 5 '14 at 12:23
  • I've edited the code to add a function (explode_rtree) that uses an rtree index, which will hopefully provide a significant performance boost. The newer function only makes one pass through the data and is more efficient with spatial tests. I did a quick test and it produced the same results as the original function, but I have not tested the code extensively. – dblanchett Aug 6 '14 at 5:46
  • That's a huge improvement - after installing Rtree (and adding "import rtree" to the script!) it took just 20 minutes to process 8700 routes producing a total of ~64000 segments to a shapefile that renders very quickly in QGIS. The output results look very plausible in terms of segment coordinates and summed totals.. and has stopped me trying to faff with a raster alternative. – Bristle6 Aug 6 '14 at 14:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.