[GIS] Shapefile to Network xml

Networkpgroutingshapefilexml

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.

Best Answer

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.

Nodes and Links

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)
Related Question