markdown Python GIS演算关连ライブラリまとめ

Posted 2021-05-06

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

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

## OSRM動作確認
```py
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へ
```py
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
```py
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を利用した地理データの可視化ライブラリ．

##### 点群表示
```py
%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
```py
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として読み込む
```py
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
```