Shapefile to Network xml

Question

I have a shapefile representing the US freight highway network available from Oak Ridge National Lab. I'm trying to prepare this as a .xml network for processing with MATSim. The output structure is very simple, something like

<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE network SYSTEM "http://www.matsim.org/files/dtd/network_v1.dtd">
<network name="VISUM export national network 2007-11-28">

<!-- ====================================================================== -->

    <nodes>
        <node id="1000" x="730065.3125" y="220415.9531" type="2" origid="1000" />
        <node id="1001" x="731010.5" y="220146.2969" type="2" origid="1001" />
        .....
    </nodes>

<!-- ====================================================================== -->

    <links capperiod="01:00:00" effectivecellsize="7.5" effectivelanewidth="3.75">
        <link id="100365" from="226" to="227" length="921.0" freespeed="33.3333333333333" capacity="5600.0" permlanes="2.0" oneway="1" modes="car" origid="183" type="10" />
    ...
    </links>
<!-- ====================================================================== -->

</network>

Most users of MATSim pull data from OpenStreetMap using osmosis, but I have to use this network for a variety of reasons.

I have seen lots of suggestions in different places as to what I might try. Some of these include:

1. v.net in GRASS
2. shp2pgsql suggested by pgRouting tutorial
3. manually with pgRouting, though this advice may be outdated?
4. writing code through the sp and rgdal libraries in R to build a node and link structure.
5. find someone with a TransCAD license to export the network ORNL provides along with the shapefile, though I don't know TransCAD exports to .xml.

It seems that with options 1 - 3, I'll then have to figure out a way to write the .xml from the database. Option 4 seems difficult, but mainly in terms of my time. I'm happy to learn something new, but I don't want to invest a day towards option 1 only to find that option 3 would have been better.

Or is there a simpler way?

---

Post-script:

Someone asked the MATSim developers if they had a tool to do the shp -> xml conversion, but the answer was that shapefiles are too wild to build one.

Answer 1

As I continued to explore yesterday, I discovered the networkx Python library, in particular its read_shp() and write_shp() functions.

import networkx
G = networkx.read_shp('linesfile.shp')
networkx.write_shp(G, './')

Got me a lines file with the original attributes and a points file with the nodes. I'm actually thrilled at the result, though there isn't a field for the node ID. Hopefully I can do this with just a spatial join.

---

The Solution

Well, I did it. Here's a reduced form of the python code I wrote. The full code with detailed comments is available in this gist

import networkx as nx
import lxml.etree as ET
G = nx.read_shp("fafnetworkLCC.shp")
G = nx.convert_node_labels_to_integers(G, first_label=0, 
        label_attribute = "coord")
# create element tree structure
network = ET.Element("network", 
    attrib={'name':"MATSim network exported from FAF shapefile."})
nodes = ET.SubElement(network, "nodes")

for i in range(len(G)):
    ET.SubElement(nodes, "node", 
            attrib={'id': str(G.nodes()[i]), 
                    'x':str(G.node[i]['coord'][0]), 
                    'y':str(G.node[i]['coord'][1]),
                    'type':"2"})

links = ET.SubElement(network, "links", 
        attrib={'capperiod': "01:00:00",
                'effectivecellsize': "7.5",
                'effectivelanewidth': "3.75"})
length = nx.get_edge_attributes(G, "MILES")
idvar  = nx.get_edge_attributes(G, "ID")

for i in range(len(G.edges())):
    startnode = G.edges()[i][0]
    endnode = G.edges()[i][1]
    ET.SubElement(links, "link", attrib={ 
        'id': str(idvar[(startnode, endnode)]),
        'from': str(startnode), 
        'to':   str(endnode), 
        'capacity': str(6000),
        'modes': "car",
        'oneway': str(1),
        'type': str(10),
        'length': str(length[(startnode, endnode)] * 1609.34)}) # convert to meters

tree = ET.ElementTree(network)

with open('network.xml', 'w') as f:
    f.write("""<?xml version="1.0" encoding="UTF-8" ?>
<!DOCTYPE network SYSTEM "http://www.matsim.org/files/dtd/network_v1.dtd">
""")
    tree.write(f, pretty_print = True)

Answer 2

Not sure if simpler, but as MATSim is written in Java, you could also write a converter in Java, if you're familiar with programming. MATSim includes a library which provides access to shape files, so it should be possible with the MATSim code to do this.

Although, looking at the actual data, it shows why MATSim probably doesn't provide a default shp-converter: There is no actual information about nodes or intersections! The polylines building the road segments just end somewhere, and one has to hope that at least another polyline starts at exactly the same coordinate in order to not have a gap in the road. Also, if two polylines intersect: do they represent a crossing, or is it a tunnel/bridge construction where vehicles have to stay on one road and cannot change to the other? (I did not check the whole dataset if this is a problem here, but I've seen such shape files where this was a problem).