Python Remind Memo - OSMnx and Folium / about Error of saving to GeoPackege type file from GeoDataFrame
Helsinki Univ. のAutomating GIS Processes 2022を参考にしながら、OSMnx(Open Street Map nx)を活用したShortest PathのSearch機能を学習しました。
その一連の学習内容をメモしておきます。
題材は京都中京区で幕末騒乱において有名になった蛤御門から菅原道真公を祀っている北野天満宮までのShortest Pathを検索し、図化するという内容です。
まずは、Google Colabでosmnx、networkx、matplotlib.pyplotという３つのlibraryをimportします。
!pip install osmnx
import osmnx
import networkx as nx
import matplotlib.pyplot as plt
次に、京都市上京区のOSMnxのgraphDataをReadします。この際、Shortest Pathについて自動車(Drive)と徒歩(Walk)を比べたいと思い立ったので、Driveに対するGraphDataとWalkに対するGraphDataをReadしたうえで、それぞれのGraphDataを図化してみます。
# Get the address of show place
PLACE_NAME = "上京区, 京都市, 京都府, 日本"
# Get networks for car
graph_car = osmnx.graph_from_place(
PLACE_NAME,
network_type='drive')
# Get networks for walk
graph_walk = osmnx.graph_from_place(
PLACE_NAME,
network_type='walk')
# Creat plot for car
osmnx.plot_graph(graph_car, edge_color='blue')
# Creat plot for walk
osmnx.plot_graph(graph_walk, edge_color='red')
次にOSMnxのgraphデータに座標系を設定します。デフォルトの座標系は緯度経度表示でUnitがDegreeになっているので、特にShortestPathAnalysisには不向きです。なので、UTM座標系に変換します。
osmnx.project_graph()メソッドを利用しますが、この際、CRSを指定しなければ、自動的に近傍のUTMのCRSをSelectしてくれます。
# Transform CRS to UTM system
# If we omit CRS parameter, OSMnx estimate CRS of neighborhood
graph_car_utm = osmnx.project_graph(graph_car)
graph_walk_utm = osmnx.project_graph(graph_walk)
# Get the nodes, edges from graph_car_utm and graph_walk_utm
nodes_car, edges_car = osmnx.graph_to_gdfs(graph_car_utm)
nodes_walk, edges_walk = osmnx.graph_to_gdfs(graph_walk_utm)
# Get the coordinates of the origin point(Hamagurigomonn)
origin = (
osmnx.geocode_to_gdf("蛤御門, 京都")
.to_crs(edges_car.crs)
.at[0, "geometry"]
.centroid # Get the centroid point
)
# Get the coordinates of destination point(Kitano Shrine)
destination = (
osmnx.geocode_to_gdf("北野天満宮, 京都")
.to_crs(edges_car.crs)
.at[0, "geometry"]
.centroid # Get the centroid point
)
print("car CRS", edges_car.crs)
print("walk CRS", edges_walk.crs)
print("origin", origin)
print("destination", destination)
Outputは以下の通りとなります。
car CRS +proj=utm +zone=53 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +type=crs
walk CRS +proj=utm +zone=53 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +type=crs
origin POINT (569296.1954866581 3875871.532623019)
destination POINT (567052.7983118674 3876677.109004879) Coodinate Reference SystemはGPSで利用されている世界測地系WGS84となっています。起点となる蛤御門、終点となる北野天満宮の座標値が表示されています。
次に、GeopandasとShaply.geometryをインポートします。
!pip install geopandas
import geopandas as gpd
import shapely.geometry
起点および終点に最も近いOSMnxのnodeを検索したあと、起点から終点までのPathのなかからShortest Pathを検索します。
# Get enach node of origin and destination in graph_car_utm
origin_node_id_car = osmnx.nearest_nodes(graph_car_utm, origin.x, origin.y)
destination_node_id_car = osmnx.nearest_nodes(graph_car_utm, destination.x, destination.y)
# Get enach node of origin and destination in graph_walk_utm
origin_node_id_walk = osmnx.nearest_nodes(graph_walk_utm, origin.x, origin.y)
destination_node_id_walk = osmnx.nearest_nodes(graph_walk_utm, destination.x, destination.y)
print("origin_node_id_car", origin_node_id_car)
print("origin_node_id_walk", origin_node_id_walk)
print("destination_node_id_car", destination_node_id_car)
print("destination_node_id_walk", destination_node_id_walk)
print("")
# Get each shortest path for car and walk
route_car = osmnx.shortest_path(graph_car_utm, origin_node_id_car, destination_node_id_car)
route_walk = osmnx.shortest_path(graph_walk_utm, origin_node_id_walk, destination_node_id_walk)
Outputは以下の通りです。
origin_node_id_car 2125132335
origin_node_id_walk 2785660326
destination_node_id_car 2376901212
destination_node_id_walk 1178978913 それぞれのNode IDが表示されます。
次がちょっと複雑だなぁと感じました。
検索したShortest PathのNodeを格納しているroute_carおよびroute_walkを参照したうえで、それぞれのNode番号が持っている座標などの情報をroute_nodes_carおよびroute_nodes_walkに格納する作業を行います。
# Get the node information(UTM coodinate, lon/lat, etc.) of the node number of shortest path
route_nodes_car = nodes_car.loc[route_car]
route_nodes_walk = nodes_walk.loc[route_walk]
route_nodes_walk
出力は以下の通りでございます。
次にそれぞれの最短経路をLineString化するとととも、GeopandasのGeoDataFrame化します。
Create a geometry for the shortest path
route_line_car = shapely.geometry.LineString(
list(route_nodes_car.geometry.values)
)
route_line_walk = shapely.geometry.LineString(
list(route_nodes_walk.geometry.values)
)
route_geom_car = gpd.GeoDataFrame(
{
"geometry": [route_line_car],
"osm_nodes": [route_car],
},
crs=edges_car.crs
)
route_geom_walk = gpd.GeoDataFrame(
{
"geometry": [route_line_walk],
"osm_nodes": [route_walk]
}
)
# Calculate the route length
route_geom_car["length_m"] = route_geom_car.length
route_geom_walk["length_m"] = route_geom_walk.length
実はこのコードにMistakeが存在します。これを解決するために２日掛かりました。
取り急ぎ、route_geom_carの内容を表示してみます。
route_geom_car
出力は以下の通りで１行だけのGeoDataFrameになっています。Route_geom_walkも座標値は異なりますが、同じような内容です。
| geometry | osm_nodes | length_m | |
|---|---|---|---|
| 0 | LINESTRING (569276.118 3875808.362, 569267.807... | [2125132335, 313317267, 433922464, 433922480, ... | 2737.902796 |
以下のコードで２つのGeoDataFrameのLineStringを図化してみると、少し形状が異なっていることが理解できると思います。
route_line_car
route_line_walk
建物を含めた背景にShortestPathを表示するために、OSMnxからBuildingInformationを読み込みます。合わせて読み込み時に座標系を統一しておきます。
buildings = osmnx.geometries_from_place(
PLACE_NAME,
{
"building": True
}
).to_crs(edges_walk.crs)
背景を作成する'contextily'モジュールと図化用の'matplotlib.pyplot'モジュールをインポートします。
!pip install contextily
import contextily
import matplotlib.pyplot as plt
ShoetestPathを図化する前に、Driveのedgeデータ（青線）とWalkのedgeデータ（赤線）を重ねて表示して、違いを確認してみました。
# Overlay edges_car and edges_walk
fig, ax = plt.subplots(figsize=(12,8))
edges_car.plot(ax=ax, linewidth=1.0, color='blue')
edges_walk.plot(ax=ax, linewidth=0.5, linestyle='--', color='red')
matplotlib.pyplotの機能を使って、CarPathとWalkPathを描き、その上に、carとWalkそれぞれのShoetestPathを太線で描き、さらに、Contextilyを使って地上物などの背景を描画します。
Contextilyは面白い機能です。知識を深めるとGISの表示に活用できるような気がします。Contectily has some kind of map, at this time I used 'contextily.providers.CartoDB.Positron' at left side graph, and 'contextily.providers.Stamen.Watercolor' at right side graph.
Next, ShortestPathのAnalysisResultをFileSaveします。
PathlibをImportしてCurrentDirectoryを検索し、さらにOsをImportしてDataSaveするDirectoryを作成します。
import pathlib
NOTEBOOK_PATH = pathlib.Path().resolve()
DATA_DIRECTORY = NOTEBOOK_PATH / "data"
import os
directory = "/content/data/"
os.makedirs(directory, exist_ok=True)
次は、CarとWalkのそれぞれのedge, route dataからInvalid Column のみを抽出した上で、Fileへのwriteのin advance processとしてAll DataをStringTypeに変換します。
# Columns with invalid values
problematic_columns_car = [
"osmid",
"lanes",
"name",
"highway",
"maxspeed",
"reversed",
"bridge",
"tunnel",
]
problematic_columns_walk = [
"osmid",
"lanes",
"name",
"highway",
"maxspeed",
"reversed",
"bridge",
"tunnel",
]
# convert selected columns to string format
edges_car[problematic_columns_car] = edges_car[problematic_columns_car].astype(str)
edges_walk[problematic_columns_walk] = edges_walk[problematic_columns_walk].astype(str)
route_geom_car["osm_nodes"] = route_geom_car["osm_nodes"].astype(str)
route_geom_walk["osm_nodes"] = route_geom_walk["osm_nodes"].astype(str)
これをそれぞれOSM_Nakagyo_car.gpkgとOSM_Nakagyo_walk,gpkgというfilenameでgeopackage typeで出力します。
# Save one layer after another
output_gpkg_car = DATA_DIRECTORY / "OSM_Nakagyo_car.gpkg"
edges_car.to_file(output_gpkg_car, layer="streets_car")
route_geom_car.to_file(output_gpkg_car, layer="route_car")
nodes_car.to_file(output_gpkg_car, layer="nodes_car")
#buildings[['geometry', 'name', 'addr:street']].to_file(output_gpkg, layer="buildings")
# Save one layer after another
output_gpkg_walk = DATA_DIRECTORY / "OSM_Nakagyo_walk.gpkg"
edges_walk.to_file(output_gpkg_walk, layer="streets_walk") # This row has unresolve error
route_geom_walk.to_file(output_gpkg_walk, layer="route_walk")
nodes_walk.to_file(output_gpkg_walk, layer="nodes_walk")
#buildings[['geometry', 'name', 'addr:street']].to_file(output_gpkg, layer="buildings")
ここでOSM_Nakagyo_walk.gpkgを作成する際にerror like belowがocureします。
>>>> ValueError: Invalid field type <class 'list'>
ChatGPT先生に尋ねたところ、どこかにlist型のDataが存在しているのではないかという回答でした。このためOriginal GeoDataFrameのDataにlisttypeのvalueがcontentしていないか、confirmしてみます。
route_geom_walk.head()
print(edges_walk.dtypes)
osmid object
oneway bool
name object
highway object
reversed object
length float64
geometry geometry
bridge object
lanes object
ref object
maxspeed object
service object
tunnel object
width object
dtype: objectroute_geom_car.head()
print(edges_walk.dtypes)
osmid object
oneway bool
name object
highway object
reversed object
length float64
geometry geometry
bridge object
lanes object
ref object
maxspeed object
service object
tunnel object
width object
dtype: object しかし、Errorのcouseであるroute_geom_walkのeach dataをviewしても、もう片方のroute_geom_carをviewしても、list type dataはexistしません。
実はここから相当な紆余曲折を経て、Errorのoriginを追求したのですが、ここではresultのみをintroduceします。
とりあえず、errorをignoreして、Folium Moduleをuseしたshow map plotにchallengeしました。（このignoreをしてとりあえずfoliumに行ってみたときにerrorを解決するヒントが見つかりました）
!pip install folium
import folium
# Read route_car data from gpkg file
map = folium.Map(location=[35.0116135, 135.7681], zoom_start=14)
# Set route_geom_car data to folium map(m)
folium.GeoJson(route_geom_car).add_to(map)
map
# initialize the folium map. Center Point is set up arround the center of Kyoto city.
map = folium.Map(location=[35.0116135, 135.7681], zoom_start=14)
# Set route_geom_car data to folium map(m)
folium.GeoJson(route_geom_walk).add_to(map)
map
当然のことながら、car routeはplotされますが、walk routeはerrorのままなのでplotはshowできません。そのときのerror massegeが次のとおりです。
>>> ValueError: Cannot transform naive geometries. Please set a crs on the object first.
geometriesを変換できないので、まず、objectのcoodinate reference systemを設定しなさいということでした。
for the confirmation, I tryed to output the coordinate reference sysytem of each GeoDataFrame data.
print(route_geom_car.crs)
print(route_geom_walk.crs)
the output is like below.
+proj=utm +zone=53 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +type=crs
NoneCoordinate reference system in route_geom_walk is 'None', but one is route_geom_car has the CRS. So I think that this condition is the course of Error.
While I write almost same code for car and walk, Why only walk side happen error. I review the code from start. I discover the difference code between car and walk.I forgot the instruction of CRS for route_geom_walk only. I try to revice the code.
# Create a geometry for the shortest path
route_line_car = shapely.geometry.LineString(
list(route_nodes_car.geometry.values)
)
route_line_walk = shapely.geometry.LineString(
list(route_nodes_walk.geometry.values)
)
route_geom_car = gpd.GeoDataFrame(
{
"geometry": [route_line_car],
"osm_nodes": [route_car],
},
crs=edges_car.crs
)
route_geom_walk = gpd.GeoDataFrame(
{
"geometry": [route_line_walk],
"osm_nodes": [route_walk]
},
crs=edges_walk.crs # <<< Add this code for improving the Error
)
# Calculate the route length
route_geom_car["length_m"] = route_geom_car.length
route_geom_walk["length_m"] = route_geom_walk.length
I try to add 'crs=edges_walk.crs'
I get good result and happenn no error without
>>> 'edges_walk.to_file(output_gpkg_walk, layer="streets_walk") # This row has unresolve error'.
I try to omit this code like below.
# Save one layer after another
output_gpkg_car = DATA_DIRECTORY / "OSM_Nakagyo_car.gpkg"
edges_car.to_file(output_gpkg_car, layer="streets_car")
route_geom_car.to_file(output_gpkg_car, layer="route_car")
nodes_car.to_file(output_gpkg_car, layer="nodes_car")
#buildings[['geometry', 'name', 'addr:street']].to_file(output_gpkg, layer="buildings")
# Save one layer after another
output_gpkg_walk = DATA_DIRECTORY / "OSM_Nakagyo_walk.gpkg"
#edges_walk.to_file(output_gpkg_walk, layer="streets_walk") # This row has unresolve error
route_geom_walk.to_file(output_gpkg_walk, layer="route_walk")
nodes_walk.to_file(output_gpkg_walk, layer="nodes_walk")
#buildings[['geometry', 'name', 'addr:street']].to_file(output_gpkg, layer="buildings")
No error has ocured.
Finally I show the plot of shortest path from Hamagurigomon to Kitano Shrine in folium from saving files.
# Read route_car data from gpkg file
route_geom_car = gpd.read_file(output_gpkg_car, layer='route_car')
route_geom_walk = gpd.read_file(output_gpkg_walk, layer='route_walk')
# initialize the folium map. Center Point is set up arround the center of Kyoto city.
map = folium.Map(location=[35.0116135, 135.7681], zoom_start=14)
# Set route_geom_car data to folium map(m)
folium.GeoJson(route_geom_car).add_to(map)
folium.GeoJson(route_geom_walk, style_function=lambda feature: {'color': 'red'}).add_to(map)
map
I want to be more beutiful map. Shortest walk route is showed by red dot thin line, and shortest car route is showed by blue straight and bold line.
# Read route_car data from gpkg file
route_geom_car = gpd.read_file(output_gpkg_car, layer='route_car')
route_geom_walk = gpd.read_file(output_gpkg_walk, layer='route_walk')
# initialize the folium map. Center Point is set up arround the center of Kyoto city.
map = folium.Map(location=[35.0116135, 135.7681], zoom_start=14)
# Set route_geom_car data to folium map(m)
folium.GeoJson(route_geom_car, style_function=lambda feature: {'color': 'blue', 'weight': 4}).add_to(map)
folium.GeoJson(route_geom_walk, style_function=lambda feature: {'color': 'red', 'weight': 2, 'dashArray': '5, 5'}).add_to(map)
map
This secion is finished here.
Thanks your cooperation!
コメント
コメントを投稿