OSRMを利用してマップマッチング

OSRM起動

$ sudo docker run -i -p 5000:5000 -v $(pwd):/data osrm/osrm-backend osrm-routed --algorithm mld  /data/kansai-latest.osrm

OSRM動作確認

def test_server():

    url = "http://" + SERVER_ADDRESS + "/route/v1/driving/13.388860,52.517037;13.397634,52.529407;13.428555,52.523219?overview=false"
    r = requests.get(url)
    result = r.json()
    
    if result['code'] != 'Ok':
        print("Error", result['code'])
        
    return result

geopandas dataframeを(lon,lat)のtupleへ

def gdf2xy_tuple_list(gdf, geom="position"):
    return [ (row[geom].coords[0][0],row[geom].coords[0][1]) for i,row in gdf.iterrows() ]

mapmatching xy tuple list -> mapmached tuple list

def request_mapmatching(xy_tuple_list):
    out_str = ""

    for point in xy_tuple_list:
        lon = point[0]
        lat = point[1]
        out_str += str(lon)+","+str(lat)+";"
        
    url = "http://" + SERVER_ADDRESS + "/match/v1/car/" + out_str[:-1]
    r = requests.get(url)
    
    result = r.json()
    
    if result['code'] != 'Ok':
        print("Error", result['code'])
        return None
        
    result_xy_tuple_list = []
    trace_points = result['tracepoints']
    for i, trace_point in enumerate(trace_points):
        if trace_point is None:
            result_xy_tuple_list.append(xy_tuple_list[i])
        else:
            loc = trace_point['location']
            result_xy_tuple_list.append((loc[0], loc[1]))
        
    return result_xy_tuple_list

Mapbox

Mapbox Mapmatching API

PythonからGETリクエストを投げる．

import requests
MAPBOX_ACCESS_TOKEN = "*"

def request_mapmatching(point_list):
    out_str = ""

    for point in point_list:
        lon = point[0]
        lat = point[1]
        out_str += str(lon)+","+str(lat)+";"
        
    url = "https://api.mapbox.com/matching/v5/mapbox/walking/" + out_str[:-1]
    url_param = {'access_token': MAPBOX_ACCESS_TOKEN }
    r = requests.get(url,params=url_param)
    
    result = r.json()
    
    if result['code'] != 'Ok':
        print("Error", result['code'])
        
    return result

Folium

Leaflet + OpenStreetMapを利用した地理データの可視化ライブラリ．

点群表示

%matplotlib inline
%config InlineBackend.figure_format='retina'
import matplotlib.pyplot as plt
import folium

m = folium.Map(location=[34.8, 135.5], zoom_start=12, tiles='cartodbpositron')
for i, row in point_gdf.iterrows():
    lon, lat = row.geometry.coords[0]
    folium.CircleMarker(location=[lat,lon],radius=2,color='#00FF00',fill=True, fill_color='#FF0000', opacity='0.9').add_to(m)

Geopandas

pandasに空間データ型を追加したライブラリ． 内部ではデータ型にShapelyを利用している．

Pandas Dataframe -> GeoPandas DataFrame

import pandas as pd
import geopandas as gpd

def convert_dataframe_to_geodataframe(df, lon_column_name, lat_column_name):
  """ Convert pandas LonLat DataFrame to GeoPandas DataFrame """
  geometry = [Point(xy) for xy in zip(df[lon_colomun_name], df[lat_colomun_name])]
  crs = {'init':'epsg:4326'}
  gdf  = gpd.GeoDataFrame(df, crs=crs, geometry=geometry)
  return gdf

PostGISのデータをGeoPanda DataFrameとして読み込む

import geopandas as gpd
import psycopg2

con = psycopg2.connect(database="", user="", password="", host="")  # set db connection

sql = """
SELECT position FROM observations
WHERE position IS NOT NULL;
"""

gdf = gpd.GeoDataFrame.from_postgis(sql, con, geom_col='positon')

geom_colで指定するカラムにNULLが存在するとエラーになるので注意．

参考：geopandas経由ではなく直接SQLを実行する

cur = con.cursor()
cur.execute(sql)

参考：pandas data frameとして読み込む

import pandas as pd
df = pd.read_sql(sql, con, index_col='id')

pyproj

座標系・測地系変換を行うライブラリ．

緯度経度（WGS84）で指定した2地点間の直線距離（meter）を求める

from pyproj import Geod
def calculate_lonlat_distance(lon1, lat1, lon2, lat2):
  q = Geod(ellps='WGS84')
  fa, ba, d = q.inv(lon1,lat1, lon2,lat2)
  return d