2017-10-04 14 views
0

C'est assez compliqué, mais pour faire court: j'ai utilisé plusieurs bibliothèques comme OSMNx pour dessiner un itinéraire entre plusieurs endroits d'une ville. Maintenant, je voudrais le convertir en un fichier shp.Shapely: tuples sans virgules avec LineString

La route est une liste complète de l'ID des nœuds. Ensuite, ces identifiants ont été utilisés pour extraire la latitude et la longitude de chaque nœud. J'ai fait tuples concaténer les coordonnées de chaque noeuds couple (un début, une arrivée) avec une boucle, comme celui-ci:

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 

Voici le résultat de l'impression (voyage) à la fin de la boucle:

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

Chaque ligne apparaît correctement. Mais quand je veux convertir dans un voyage LineString de bien fait ... Et il retourne ceci:

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) 

Par conséquent, je ne peux pas le convertir dans un shp en utilisant 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

Je ne comprends pas pourquoi les tuples sont convertis sans une virgule entre la latitude et la longitude. De plus, il y a plusieurs duplications (les deuxièmes coordonnées de la troisième ligne sont les premières coordonnées de la quatrième ligne, etc ...) et peut-être cela pourrait être une source d'erreur pour le futur shp.

Merci d'avance!

+0

Ce que vous voyez quand vous faites 'print (final_journey)' est la représentation de [Texte bien connu] (https://en.wikipedia.org/wiki/Well-known_text) de votre ligne. Il n'y a rien de mal à cela (comme 'tuples sans virgule'), c'est juste comment bien afficher les géométries dans l'interpréteur. – mgc

Répondre

1

Je ne pense pas qu'obtenir des coordonnées de nœuds, puis les reliant ensemble est la meilleure chose que vous pouvez faire. Et si la rue n'est pas droite? OSMnx vous donne la géométrie exacte des rues. pour extraire la géométrie des nœuds, la meilleure solution est expliquée here. Mais vous avez besoin de la géométrie des rues. Comme il peut y avoir plus d'une arête entre deux nœuds, ce n'est pas toujours simple à faire. Je crois que ox.get_route_edge_attributes() devrait être capable de faire cela, et en effet cela fonctionne très bien si vous demandez un autre attribut (par exemple highway), mais pas pour extraire geometry des bords. La raison (je suppose) est que tous les bords de G n'ont pas geometry, mais si vous obtenez gdf_edges du réseau, vous avez toujours la géométrie de chaque bord. Ce qui suit est le travail autour j'ai trouvé:

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

Pour obtenir le GeoDataFrame des noeuds dans la route:

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

enter image description here

et pour obtenir la géométrie des bords d'abord définir le tuple u, v valeurs comme index du gdf_edges, puis loc pour sélectionner le GeoDataFrame du chemin:

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

vous pouvez l'enregistrer sur shapefile:

output_edges.to_file() 

Une remarque importante: comme je l'ai déjà dit, il peut y avoir plus d'un bords entre deux nœuds.Cela signifie que pour identifier une arête de manière unique, u et v (début et fin du bord) ne suffisent pas, vous avez également besoin de key, qui est automatiquement ajouté par OSMnx et vous le trouvez à la fois dans graph et gdf_edges pour chaque arête. Donc, si vous utilisez le code ci-dessus, sachez qu'il vous donnera all les bords (s'il y en a plus d'un) entre les nœuds. une vérification rapide est: len(np.nonzero(output_edges['key'])[0]). s'il n'y a pas de bords parallèles, ce sera sûrement zéro. sinon, cela signifie qu'il y a des arêtes parallèles.

MISE À JOUR

OSMnx a une fonction pour enregistrer shapefile de GeoDataFrame:

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

ajoutant schéma à to_file() semble fonctionne aussi:

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) 
+0

Salut Alireza! Juste testé votre solution et cela a fonctionné bien jusqu'à l'étape output_edges.to_file() (je ne suis pas sûr que c'est la bonne variable, je pensais que c'était mieux avec output_lines). Dois-je soumettre un schéma comme celui de ma question? Merci d'avance ! – Raphadasilva

+0

@Raphadasilva consultez la mise à jour. –

+0

Merveilleux, je valide ta réponse! Merci encore – Raphadasilva