0
votes

it's quite complicated, but long story short : I've used several libraries like OSMNx to draw a route between several spots of a city. Now I would convert it as a shp file.

The route is a list full of nodes' id. Then these id were used to extract latitude and longitude of each node. I made tuples to concatenate the coordinates of each nodes couple (one start, one arrival) with a for loop, like this :

journey = [] 
# previous list will contain tuples with coordinates of each node

for node1, node2 in zip(route[:-1], route[1:]):
    parcours.append(tuple((G.node[noeud1]['x'], G.node[noeud1]['y']))) # we create a tuple with coordinates of start's node
    parcours.append(tuple((G.node[noeud2]['x'], G.node[noeud2]['y']))) # then we make the same for the arrival node

Here's the result of print(journey) at the end of the loop :

[(6.15815, 48.6996136), (6.1629696, 48.7007431), (6.1629696, 48.7007431), [...], (6.1994411, 48.6768434), (6.1994411, 48.6768434), (6.1995322, 48.6767583)]

Each tuple appears correctly. But when I want to convert journey in a shapely's LineString... And it returns this :

from shapely.geometry import LineString
final_journey = LineString(journey)
print(final_journey)
LINESTRING (6.15815 48.6996136, 6.1629696 48.7007431, 6.1629696 48.7007431, 6.1630717 48.7002871, [...], 6.1991794 48.677085, 6.1994411 48.6768434, 6.1994411 48.6768434, 6.1995322 48.6767583)

Consequently, I can't convert it in a shp using fiona :

import fiona
schema = {
    'geometry': 'Polygon',
    "properties": {'id': 123}
}

with fiona.open('test.shp', 'w', 'ESRI Shapefile', schema) as c:
    c.write({
        'geometry': mapping(trace)
    })

--------------------------------------------------------------------------- TypeError Traceback (most recent call last) in () 4 } 5 ----> 6 with fiona.open('test.shp', 'w', 'ESRI Shapefile', schema) as c: 7 c.write({ 8 'geometry': mapping(trace)

/usr/local/lib/python3.5/dist-packages/fiona/init.py in open(path, mode, driver, schema, crs, encoding, layer, vfs, enabled_drivers, crs_wkt) 173 c = Collection(path, mode, crs=crs, driver=driver, schema=this_schema, 174 encoding=encoding, layer=layer, vsi=vsi, archive=archive, --> 175 enabled_drivers=enabled_drivers, crs_wkt=crs_wkt) 176 else: 177 raise ValueError(

/usr/local/lib/python3.5/dist-packages/fiona/collection.py in init(self, path, mode, driver, schema, crs, encoding, layer, vsi, archive, enabled_drivers, crs_wkt, **kwargs) 154 elif self.mode in ('a', 'w'): 155 self.session = WritingSession() --> 156 self.session.start(self, **kwargs) 157 except IOError: 158 self.session = None

fiona/ogrext.pyx in fiona.ogrext.WritingSession.start (fiona/ogrext2.c:16207)()

TypeError: argument of type 'int' is not iterable

I don't understand why tuples are converted without a comma between latitude and longitude. In addition, there are several duplication (the second coordinates of the third line is the first coordinates of the fourth line, etc...) and maybe it could be source of error for the future shp.

Thanks in advance !

1
What you are seeing when you do print(final_journey) is the Well Known Text representation of your line. There is nothing wrong with this (such as 'tuples without comma'), it's just how shapely display the geometries in the interpreter.mgc

1 Answers

1
votes

I don't think getting coordinates of nodes and then connecting them together is the best thing you can do. what if the street is not straight? OSMnx gives you the exact geometry of streets. to extract geometry of nodes, the better solution is explained here. But you need geometry of streets. Since there might be more than one edge between two nodes, this is not always simple to do. I believe ox.get_route_edge_attributes() should be able to do that, and indeed it works great if you ask for another attribute (e.g. highway), but not for extracting geometry of edges. The reason (I guess) is that not all edges in G have geometry, but if you get gdf_edges of the network, then you always have the geometry of each edge. The following is the work around I found:

gdf_nodes, gdf_edges = ox.graph_to_gdfs(G)
path = nx.shortest_path(G, G.nodes()[0], G.nodes()[1])

To get the GeoDataFrame of nodes in the route:

output_points = gdf_nodes.loc[path]
output_points.plot(color='red')

enter image description here

and to get the geometry of edges first set the tuple of u,v values as index of the gdf_edges, then loc to select the GeoDataFrame of the path:

gdf_edges.index = gdf_edges.apply(lambda row: (row.u, row.v), axis=1)
output_lines = gdf_edges.loc[list(zip(path[:-1], path[1:]))]
output_lines.plot(color='red')

enter image description here

you can then save it to shapefile:

output_edges.to_file()

One important remark: as I already said, there can be more than one edges between two nodes. This means that to uniquely identify an edge, u and v (start and end of the edge) are not enough, you also need key, which is added by OSMnx automatically and you find it in both graph and gdf_edges for every edge. so if you use the above code, be aware that it will give you all the edges (if there are more than one) between nodes. a quick check is: len(np.nonzero(output_edges['key'])[0]). if there is no parallel edges, this will surely be zero. if not, it means that there are parallel edges.

UPDATE

OSMnx has a function to save shapefile from GeoDataFrame:

p = '/path/to/destination/'
ox.save_gdf_shapefile(output_lines, 'output', p)

adding schema to to_file() seems to works too:

sch = {'geometry': 'LineString',
       'properties': {'key': 'int',
                        'u': 'int',
                        'v': 'int'}}

test = gpd.GeoDataFrame(output_lines[['u', 'v', 'key']], 
                        geometry=output_lines['geometry'], crs='+init=EPSG:3740')

test.to_file('/path/to/destination/test.shp', schema=sch)