3

I want to export a directed weighted graph from a shapefile. I used read_shp function of the Networkx package to export the directed graph which perfectly matches my needs. But I am unable to calculate the length of each edge as line geometries are simplified into start and end coordinates in the output of Networkx. How can I calculate the exact length of the edges (linestrings)?

  • Help says that attributes are preserved, so just try something like G[f][t]['length'] if there is a field called "Length" in the shapefiles table – FelixIP Oct 7 '16 at 2:49
3

see Inaccurate output (missing features) while reading a shapefile into networkx) and nx_spatial does not read all shapefile features

Networkx generate a networkx.DiGraph with nodes without duplicates.

import networkx as nx
G = nx.read_shp('edges_length_stac.shp'

The original LineStrings and the resulting nodes of the graph

enter image description here

And the calculated distance is always between the blue nodes.

But the original geometry is still present in the edge data

First edge

 first = G.edges(data=True)[0]
 print first
 ((203153.19849954147, 89071.43217150889), (202791.36231525266, 89985.30048311011), {'ShpName': 'edgse_length', 'Json': '{ "type": "LineString", "coordinates": [ [ 203153.198499541467754, 89071.432171508888132 ], [ 202833.112644209060818, 89521.407939150099992 ], [ 202684.667030141863506, 89781.187763767709839 ], [ 202791.362315252656117, 89985.300483110113419 ] ] }', 'weight': 982.8940508567346, 'Wkb': '\x00\x00\x00\x00\x02\x00\x00\x00\x04A\x08\xcc\x89\x96\x86\xedw@\xf5\xbe\xf6\xea,\xac\x0fA\x08\xc2\x88\xe6\xb2\x01\xd0@\xf5\xdb\x16\x86\xeb3\xc7A\x08\xbd\xe5V\x13\xe6&@\xf5\xebS\x01\x14\x94\x9bA\x08\xc1:\xe6\x05\x8a\x08@\xf5\xf8\x14\xce\xc7`\xaf', 'id': None, 'Wkt': 'LINESTRING (203153.198499541467754 89071.432171508888132,202833.112644209060818 89521.407939150099992,202684.667030141863506 89781.187763767709839,202791.362315252656117 89985.300483110113419)'}) 

The two first elements of the resulting list are the nodes coordinates and the distance between is

distance = math.sqrt(sum([(a - b) ** 2 for a, b in zip(first[0],first[1])]))
print distance
982.8940508567346

But what interests us is the third element which is a dictionary.

print first[2].keys()
['ShpName', 'Json', 'weight', 'Wkb', 'id', 'Wkt']
# in GeoJSON format
print first[2]['Json']
{ "type": "LineString", "coordinates": [ [ 203153.198499541467754, 89071.432171508888132 ], [ 202833.112644209060818, 89521.407939150099992 ], [ 202684.667030141863506, 89781.187763767709839 ], [ 202791.362315252656117, 89985.300483110113419 ] ] }
# in WKT format
print first[2]['Wkt']
LINESTRING (203153.198499541467754 89071.432171508888132,202833.112644209060818 89521.407939150099992,202684.667030141863506 89781.187763767709839,202791.362315252656117 89985.300483110113419)

And we can calculate the original length with Shapely for example.

from shapely.geometry import shape
import json
print shape(json.loads(first[2]['Json'])).length
1081.7261468150266

Now more simply

for i in list(G.edges_iter(data='Json')):
    print i
((203153.19849954147, 89071.43217150889), (202791.36231525266, 89985.30048311011), '{ "type": "LineString", "coordinates": [ [ 203153.198499541467754, 89071.432171508888132 ], [ 202833.112644209060818, 89521.407939150099992 ], [ 202684.667030141863506, 89781.187763767709839 ], [ 202791.362315252656117, 89985.300483110113419 ] ] }')
((203992.8440041091, 89256.9891890929), (204192.3177980119, 90027.0508120665), '{ "type": "LineString", "coordinates": [ [ 203992.844004109094385, 89256.989189092899323 ], [ 204215.512425209890353, 89479.657610193695291 ], [ 204080.983587461494608, 89785.82668920730066 ], [ 204192.317798011907144, 90027.050812066503568 ] ] }')
((204011.3997058675, 88950.8201100793), (203403.70047327987, 90050.24543926452), '{ "type": "LineString", "coordinates": [ [ 204011.399705867486773, 88950.820110079293954 ], [ 203487.201131192676257, 89270.90596541170089 ], [ 203417.617249598668423, 89669.853553217297303 ], [ 203682.035999655898195, 89883.244123438911629 ], [ 203403.700473279866856, 90050.24543926451588 ] ] }')

New

If you want a Planar Graph in networkx you should not use nx.read_shp()

The shapefile and the nodes with nx.read_shp()

enter image description hereenter image description here

# combine the lines of the shapefile
import fiona
from shapely.geometry import shape
lines =[shape(line['geometry']) for line in fiona.open("test.shp")]

Now eliminate the intersections (a graph is planar if it can be drawn in a plane without graph edges crossing)

from shapely.ops import unary_union
result = unary_union(lines)

Use the segments of the resulting lines

G = nx.Graph()
import itertools
for line in result:
   for seg_start, seg_end in itertools.izip(list(line.coords),list(line.coords)[1:]):
       G.add_edge(seg_start, seg_end)

Resulting nodes

enter image description here

And now, you have a complete planar graph with the real distances between the nodes but without the faces

faces = nx.cycle_basis(G)
for face in faces:
    print Polygon(face)
POLYGON ((165871.1317836872 116667.8610861057, 166126.0509806138 116414.9028060785, 166514.6755202576 116900.6834806333, 166324.6383829749 116966.4195973035, 166031.229247432 117078.2517678177, 165871.1317836872 116667.8610861057))

enter image description here

  • Thanks for your answer. But I observed your suggested answers for other mentioned questions. Do you think that the NetworkX is not suitable for extracting road network form shapefile (as the NetworkX module do not support Planar Graph)? – user302787 Oct 8 '16 at 20:05
  • look above in New – gene Oct 9 '16 at 16:15

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.