交通情報学特論 第7回「PostgreSQL + PostGIS + QGIS による公共交通データ分析 3」講師：伊藤昌毅

東京大学 大学院情報理工学系研究科 創造情報学専攻 「交通情報学特論」 第7回 2023年5月31日 弥生キャンパス I-REF棟 Hilobby PostgreSQL+PostGIS+QGISによる 交通データ分析 (3) Join応用、ウィンドウ関数 東京大学 大学院情報理工学系研究科 附属ソーシャルICT研究センター 創造情報学専攻兼担 伊藤昌毅

感想・コメントありがとうございました • 出席確認を兼ねているので、履修者は必ず提出してください – 次回の水曜日23:59:59を締切としますが、締切後も提出できるはず

PostGISとjoinの応用

QGIS上のデータをPostgreSQLに持っていく • 第2回授業で実施した国勢 調査人口データを利用 • Qiita記事の後半「PostGIS に読み込ませるためのデー タ書きだし」でも解説 – https://qiita.com/niyalist/item s/d70f471c259211aa1554

• https://qiita.com/niyalist/items/d70f471c259211aa1554

授業内ではバックアップデータから復元 • x https://www.dropbox.com/s/y7p3xg1aybyrlym/population.backup?dl=0

• https://www.dropbox.com/s/y7p3xg1aybyrlym/population.backup?dl=0

国勢調査データをPostgreSQL dump形式 で保存 • Schema を base, create schema をNO に • ファイル名がtable名になる

コマンドラインからデータインポート • それぞれ書き出したファイルを置いたフォルダをカレントディ レクトリに設定したうえで • Windows C:¥census-data¥population>"C:¥Program Files¥PostgreSQL¥15¥bin¥psql.exe" -U postgres -d TranspoCensus -f population.sql • Mac $ psql -d TranspoCensus -f population.sql – 必要に応じてユーザ名などを指定

駅勢圏人口を算出 • 駅から半径800mの円を求め、その中の人口を算出 – 駅勢圏の距離については都市か地方かなどによって異なる値が用いられる • GISにおけるBuffering – ある地物に対して、指定の距離の領域の図形を作成 – 点ならば円、線ならば幅のある線など

800mをどう指定するか？ • GISは一般的に緯度経度データ を扱う – →縦横方向で1の長さが異なり、緯度 によっても異なるのでこの座標系のま までは正円は扱いづらい • メートル座標系を使いたい場合 は、UTM図法で投影してデー タ処理するのがよい – WGS84 UTM 54N： EPSG:32654

GISでよく使う空間座標系 その1 • 緯度経度による座標系（単位： 度） – 緯度経度の数字だけでは1点を示せない。必ず準拠し ている測地系とともに示す必要がある – 測地系： 地球の形状を定義したもの – 緯度経度のままでは、長さや面積の計算は難しい • 緯度・経度それぞれで1度の長さが異なる • 高緯度になると経度1度が短くなる • 例： – WGS84： GPSで用いられる測地系 EPSG:4326 – JGD2000： 日本測地系2000 EPSG:4612 • 実用上はWGS84と共通（世界測地系とも呼ばれる） • 東日本大震災による地殻変動を考慮したJDG2011も用いられ る – 旧日本測地系： EPSG:4301 • 明治時代の測量成果に基づく。古いシステムやデータはこれに 基づく場合がある https://www.gsi.go.jp/KIDS/PAMPHLET/p9.htm

• https://www.gsi.go.jp/KIDS/PAMPHLET/p9.htm

GISでよく使う空間座標系 その2 • 投影座標系（単位： メートルなど） – 地図を平面に表現したり、データ処理を容易にするため、緯度経度を2次元に投影した 座標系 – 距離や面積の計算が容易だが、投影による歪みが存在することに注意。ある投影法は、 限られた領域において有効なことが多い • 例： – Pseudo-Mercator: Google MapsなどWeb地図に使われる投影法 EPSG:3857 – UTM（ユニバーサル横メルカトル）図法 • 経度6度毎に縦に分割しその範囲を平面に投影 • 東京周辺なら54帯（EPSG:32654） – 平面直角座標系 • 日本国内の測量などに用いられる座標系 • 都道府県などにもとづき日本全体を19分割 日本周辺のUTMゾーン https://ja.wikipedia.org/wiki/ユニバーサル横メルカトル図法

• https://ja.wikipedia.org/wiki/

PostGISによるBuffering処理 • ST_Transformで一度 EPSG:32654に変換し、ST_Bufferで 800mの円を作成、更にそれをST_TransformでWGS84に戻す Select *, ST_Transform(ST_Buffer(ST_Transform(center, 32654), 800), 4326) as buffer from analysis.station_master where city is not null ; SQL-7-1 • city is not null を指定することで、市町村境界データとのjoin が成功した駅のみを抽出

人口を面積で按分 • 5次メッシュ（250mメッシュ）の人口を面積で按分 駅勢圏の形状 重複部分を持つ5次メッシュ 重複部分を切り出したメッシュ

各領域に対して、国勢調査メッシュをjoin select * from( select *, ST_Transform(ST_Buffer(ST_Transform(center, 32654), 800), 4326) as buffer from analysis.station_master where city is not null SQL-7-1 )as data1 inner join base.population as pop on ST_Intersects(data1.buffer, pop.wkb_geometry) ; SQL-7-2 • ST_intersects: 2つの図形に共通部分があればTrue

例：あざみ野駅は45行の5次メッシュ区域とマッチ 5次メッシュのエリア （約250m四方） 駅名 （同じものが45行） 半径800mの円 （同じものが45行） メッシュ内人口 45行続く

• x

人口を面積で按分する • 1. メッシュの面積を求める … α – ST_Area(図形) – 方法1: 投影座標系に事前に変換する – 方法2： geom::geography としてgeography 型にキャストすることで平方 メートル単位の面積を求められる • 2. bufferと重なる部分を抽出 – ST_Intersection（図形1, 図形2） • 3. 重なった部分の面積を求める … β • 4. 比率を求め人口と掛ける – population × β/α

例： あざみ野駅エリアの按分した人口 • x select name, population, ST_Area(wkb_geometry::geography) as area, ST_Area(ST_Intersection(buffer, wkb_geometry)::geography) as intersection_area, ST_Area(ST_Intersection(buffer, wkb_geometry)::geography)/ST_Area(wkb_geometry::geography) as rate, population * ST_Area(ST_Intersection(buffer, wkb_geometry)::geography) / ST_Area(wkb_geometry::geography) as partial_population, ST_Intersection(buffer, wkb_geometry) from 略 where name='あざみ野' ; 参考： 250m*250m = 62,500㎡ SQL-7-2改1

駅単位で集計 • 駅の primary key である駅名+県名+市町村名で group by • 人口データは、partial populationの合計を算出

select name, prefecture, city, sum(partial_population) as population from( select name, prefecture, city, population * ST_Area(ST_Intersection(buffer, wkb_geometry)::geography) / ST_Area(wkb_geometry::geography) as partial_population from( select *, ST_Transform(ST_Buffer(ST_Transform(center, 32654), 800), 4326) as buffer from analysis.station_master where city is not null )as data1 inner join base.population as pop on ST_Intersects(data1.buffer, pop.wkb_geometry) SQL-7-1 SQL-7-2改-2 )as data2 group by name, prefecture, city ; SQL-7-3

仕上げ：テーブルに保存 • 集計後にデータとして残すために、center, geom, buffer でも group by に含む – 元々の analysis.station_master を join すればいいのだけど… • create table でテーブルに出力

drop table if exists station_population; create table station_population as select name, prefecture, city, sum(partial_population) as population, center, geom, buffer from( select name, prefecture, city, population * ST_Area(ST_Intersection(buffer, wkb_geometry)::geography) / ST_Area(wkb_geometry::geography) as partial_population, center, geom, buffer from( select *, ST_Transform(ST_Buffer(ST_Transform(center, 32654), 800), 4326) as buffer from analysis.station_master where city is not null )as data1 inner join base.population as pop on ST_Intersects(data1.buffer, pop.wkb_geometry) )as data2 group by name, prefecture, city, center, geom, buffer ; SQL-7-3改

ボロノイ分割の併用 • ボロノイ分割：その位置から 最も近い点によってエリアを 分割 • ChatGPTに聞いてみた→惜 しいけど考え方は参考になる

まずはボロノイ分割 select (ST_Dump(ST_VoronoiPolygons(ST_Collect(center)))).geom AS voronoi from analysis.station_master where city is not null • 実際に用いるときは以下の関数の働きも理解すること – ST_Collect: 集約関数として全ての要素を集計 – ST_Dump: 出力された単一のgeometry collectionを個々の要素に分解 SQL-7-4

UTM図法に投影してからボロノイ分割 select (ST_Dump (ST_Transform( ST_VoronoiPolygons( ST_Collect( ST_Transform(center, 32654))), 4326))).geom AS voronoi from analysis.station_master where city is not null SQL-7-4改 • ST_Transform でUTM54系に変換→ボロノイ分割→WGS84に 戻す

• x

create table analysis.station_voronoi_buffer as with voronoi as( select (ST_Dump(ST_Transform(ST_VoronoiPolygons(ST_Collect(ST_Transform(center, (ST_Dump(ST_VoronoiPolygons(ST_Collect(center)))).geom AS voronoi 32654))), from 4326))).geom AS voronoi from analysis.station_master where analysis.station_master where city is not null ), city is not null buffer as ( ), buffer as ( select select *, ST_Transform(ST_Buffer(ST_Transform(center, 32654), 800), 4326) as buffer *, from ST_Transform(ST_Buffer(ST_Transform(center, 32654), 800), 4326) as buffer from analysis.station_master where analysis.station_master where city is not null ) city is not null )select select name, prefecture, name, city, prefecture, center, city, geom, center, ST_Intersection(voronoi, buffer) as buffer geom, ST_Intersection(voronoi, buffer) as buffer from buffer as b inner join voronoi as v on ST_Contains(v.voronoi, b.center); b.center) SQL-7-4改 SQL-7-1

修正版

修正版

修正版

動的データの収集・分析 （Window関数）

GBFS • シェアモビリティの国際標準 フォーマット – シェアサイクル、シェアキック ボードなどのポートデータが公開 されている • 日本では2022年よりドコモ バイクシェア、ハローサイク リングのデータが公開 https://www.odpt.org

• https://www.odpt.org/

主なデータ項目: それぞれjsonファイルに収納 vehicle_types.json station_information.json, • 提供されている車両の種 station_status.json 類 • ポート情報 – 自転車、キックボードなどの 種類 – 動力の種類やバッテリ残量 – 緯度経度、現在の車両台数、 返却出来る台数など – 乗り捨て自由の場合は free_bike_statusとして車両 の緯度経度などを提供 station_regions.json • サービス提供エリア system_pricing_plans.json • 利用料金 • サービスの利用者向け情報提供が目的 • シェアリングサービスの現在の状況を表現 （過去の情報は残らない）

https://api-public.odpt.org/api/v4/gbfs/hellocycling/station_information.json
（大量のポート情報のひとつを抜粋）
{
"lat": 35.707252,
ポート名、緯度経度、住所
"lon": 139.777587,
"name": "新御徒町ステーション",
"address": "東京都台東区台東４丁目３１－１",
"station_id": "17",
"rental_uris": {
"ios": "https://www.hellocycling.jp/app/port/detail/17?referrer=odpt",
"web": "https://www.hellocycling.jp/app/port/detail/17?referrer=odpt",
"android": "https://www.hellocycling.jp/app/port/detail/17?referrer=odpt"
},
"parking_hoop": false,
利用する際のURL
"parking_type": "street_parking",
"contact_phone": "+8105038218282",
"is_charging_station": false,
"vehicle_type_capacity": {
"num_bikes_now": 0,
"num_bikes_limit": 8,
"num_bikes_parkable": 8,
利用可能/返却可能な
"num_bikes_rentalable": 0
自転車の数
}
},

• https://api-public.odpt.org/api/v4/gbfs/hellocycling/station_information.json

ChatGPTを使って自動ダウンロードプログラム 作成方法を検討 • とりあえず聞くと大枠を教えて くれる – URLを直せばとりあえず使える • 更に修正するとしたら… – ダウンロードしたJSON形式をCSVに変 換してから保存 • Pandasを利用 – ファイル名にタイムスタンプを入れて、 重複がなくなるように – ダウンロードの度にデータを追記

import os import time import requests import pandas as pd from datetime import datetime service_name='HelloCycling' station_info_url = 'https://api-public.odpt.org/api/v4/gbfs/hellocycling/station_information.json' station_status_url = 'https://api-public.odpt.org/api/v4/gbfs/hellocycling/station_status.json' sleep_time = 290 def download_data(url): response = requests.get(url) data = response.json() return data def station_info_to_csv(data, filename): # "stations"キーの中身をDataFrameに変換 df = pd.DataFrame(data['data']['stations']) # タイムスタンプ列を追加（UNIXタイムスタンプを日時に変換） timestamp = datetime.fromtimestamp(data['last_updated']).strftime('%Y-%m-%d %H:%M:%S’) df['timestamp'] = timestamp now = datetime.now().strftime('%Y-%m-%d %H:%M:%S’) df['download_time'] = now df['service']=service_name df.to_csv(filename, index=False) def station_status_to_csv(data, filename): # "stations"キーの中身をDataFrameに変換 df = pd.DataFrame(data['data']['stations']) # vehicle_types_availableとvehicle_docks_availableを展開 df['vehicle_types_available_count'] = df['vehicle_types_available'].apply(lambda x: x[0]['count'] if x else None) df['vehicle_types_available_vehicle_type_id'] = df['vehicle_types_available'].apply(lambda x: x[0]['vehicle_type_id'] if x else None) df['vehicle_docks_available_count'] = df['vehicle_docks_available'].apply(lambda x: x[0]['count'] if x else None) df['vehicle_docks_available_vehicle_type_ids'] = df['vehicle_docks_available'].apply(lambda x: ', '.join(x[0]['vehicle_type_ids']) if x else None)

# 不要な列を削除
df = df.drop(['vehicle_types_available', 'vehicle_docks_available'], axis=1)
# UNIXタイムスタンプを日時に変換
df['last_reported'] = pd.to_datetime(df['last_reported'], unit='s’)
# 時間をUTCからJSTに変換
df['last_reported'] = df['last_reported'].dt.tz_localize('UTC').dt.tz_convert('Asia/Tokyo’)
# タイムゾーン情報を削除
df['last_reported'] = df['last_reported'].dt.tz_localize(None)
# タイムスタンプ列を追加（UNIXタイムスタンプを日時に変換）
timestamp = datetime.fromtimestamp(data['last_updated']).strftime('%Y-%m-%d %H:%M:%S’)
df['timestamp'] = timestamp
now = datetime.now().strftime('%Y-%m-%d %H:%M:%S’)
df['download_time'] = now
df['service']=service_name
df.to_csv(filename, mode='a', header=not os.path.exists(filename), index=False)

#最初に station_infomationを1回ダウンロード
json_data = download_data(station_info_url)
# JSONデータをCSVに変換して保存（ファイル名に日付含む）
now = datetime.now().strftime('%Y-%m-%d-%H-%M-%S')
station_info_to_csv(json_data, f'station_information_{now}.csv')
csv_file_name = f'station_status_{now}.csv'
#その後、 station_status を約5分間隔でダウンロード
while True:
# URLからJSONデータをダウンロード
json_data = download_data(station_status_url)

# JSONデータをCSVに変換して保存（ファイル名に現在の日時を含む）
station_status_to_csv(json_data, csv_file_name)
# 5分(300秒)待つ
time.sleep(sleep_time)

クローラーを作るならば • 無限ループするプログラムを作成(今回の手法） – 利点： 簡単に試せる。実行間隔が比較的安定する – 欠点： PC再起動時、トラブル時の対応が必要、PCを起動し続ける必要 • OSの定期タスク呼び出し機能を利用 – Linux, macOS, UNIX: cron – macOS: launchd – Windows: タスクスケジューラ • クラウドを利用 – AWS: Lambda + EventBridge

シェアモビリティデータをバックアップデータから 復元 • x https://www.dropbox.com/s/z20p7ovlr36vwtb/sharemobility05-25-26.backup?dl=0

• https://www.dropbox.com/s/z20p7ovlr36vwtb/sharemobility-05-25-26.backup?dl=0

ハローサイクリング・ドコモバイクシェアのポート位置

演習用データ • ハローサイクリング、ドコモバイクシェアのデータのクロール 結果 • station_information.json – クロール開始時に一度取得。ポート位置のマスターデータとして • station_status.json – クロール開始時から定期定期に取得 – 周期: 290秒（更新が5分間隔なので漏らさないため） • それぞれのデータをCSV化しPostgreSQLに読み込ませる

例：文京総合体育館 select • x * from sharemobility.station_information_docomo where station_id = '00010244' select * from sharemobility.station_status_docomo where station_id = '00010244' order by last_reported

データを眺めてみる • station_information – ポートの名称、位置、容量など – マスターテーブル • station_status – 約5分おきに利用可能な自転車数、返却可能な自 転車数が並ぶ • 合計は常に15（ポートのキャパシティ） – トランザクションテーブル

貸出可能な自転車数の変化 18 16 14 • x 12 10 8 6 4 2 0 5/25/23 0:00 5/25/23 4:48 5/25/23 9:36 5/25/23 14:24 5/25/23 19:12 5/26/23 0:00 5/26/23 4:48 5/26/23 9:36 5/26/23 14:24 5/26/23 19:12 5/27/23 0:00

利用数、返却数の推定 • 5分おきの自転車台数を前後で比較し判定 – 増：自転車が返却された – 減：自転車が利用された – （返却と利用が同時に起こることは今回は想定しない） 5/25 8:14頃 1台利用された 5/25 8:19頃 1台返却された 5/25 8:38頃 1台利用された

SQLで考える→ウィンドウ関数 • select 節の関数 – 同じ行の中での演算のみ・前後の行は見ていない • 集約関数 – group by などの条件に従って1行にまとめる – 今回は1行にまとめたいわけではない • ウィンドウ関数（分析関数） – ある範囲の行全体を見たうえで値を計算する関数 • ある範囲内で何番目に大きいか？ • 前後の行の値は？その差は？ – 集約関数に似ているが、複数行を1行にまとめない

データ処理の方針 • • • • • • 同じstation_id のデータを last_repoted 順に並べる 1行前の num_bikes_available を取得 現在の num_bikes_available と差を計算 負の場合: 自転車が利用された 正の場合： 自転車が返却された ある期間で集計

ウィンドウ関数 select station_id, last_reported, num_bikes_available, lag(num_bikes_available, 1) over(partition by station_id order by last_reported) as pre_bikes_avilable from sharemobility.station_status_docomo where station_id = '00010244' order by last_reported • over() の中にウィンドウ関数が関知する範囲について記述 – 今回は、station_id 単位で区切ることと、 last_repoted の順番を考慮することを指定

実行結果 • 1行前の num_bikes_available を得ている • 各 station_id の先頭 行においては、nullに なる 先頭はnull

ウィンドウ関数の例 関数名 count sum 意味 グループ内の行数を数える グループ内の合計や累積和をとる avg rank row_number グループ内の平均をとる グループ内をソートして順位を付ける（重複あり） グループ内をソートして順位を付ける（重複無し） lag lead グループ内をソートして前の行の値を取る グループ内をソートして後の行の値を取る 一部は集約関数と同じであることに注意 ※ 青木峰郎 10年戦えるデータ分析入門 p.170より

実際には一気に bike_diff を計算する select num_bikes_available - lag(num_bikes_available, 1) over(partition by station_id order by last_reported) as bike_diff, * from sharemobility.station_status_docomo where station_id = '00010244' order by last_reported このときも先頭はnull

正負に応じて返却、利用をカウントする select • x case when bike_diff < 0 then -bike_diff else 0 end as bike_used, case when bike_diff > 0 then bike_diff else 0 end as bike_returned, * from( select num_bikes_available - lag(num_bikes_available, 1) over(partition by station_id order by last_reported) as bike_diff, * from sharemobility.station_status_docomo )as data1

station_id ごとに集計 select station_id, sum(bike_used) as bike_used, sum(bike_returned) as bike_returned, min(last_reported) as data_start, max(last_reported) as data_end, count(*) as data_num from( select case when bike_diff < 0 then -bike_diff else 0 end as bike_used, case when bike_diff > 0 then bike_diff else 0 end as bike_returned, * from ( select num_bikes_available - lag(num_bikes_available, 1) over(partition by station_id order by last_reported) as bike_diff, * from sharemobility.station_status_docomo )as data1 ) as data2 group by station_id

マスターテーブルをjoin select * from( 前掲SQL ) as data3 inner join sharemobility.station_information_docomo as info on data3.station_id = info.station_id • create table するためには、重複の列を削除する必要 • inner join によって消えてしまう行が数行あるが今は気にしない – おそらくポートの新設などで、マスターデータのタイミングでは存在しなかったデータを収集している

完成したテーブルをDBに保存

SQL: ここまでで分析に必要な項目はある程度網羅 • 第3回 – 入門（select … from）、whereによる絞り込み、group by による集計 • 第4回 – union all による表の結合、サブクエリ（入れ子） – create table による出力から新しいテーブル作成 – PostGISによる空間データ処理 • クラスタリング、包含関係判定 – QGISとの連携 • 第6回 – Join 単純なjoinからPostGISを用いたJoin – データベースの正規化（テーブル分割） • 第7回 – 実用的なJoin – ウィンドウ関数 – with 節

更に学ぶために • Ｘ https://zenn.dev/boiledorange73/books/caea8d4c77dbba2e23a0

• https://zenn.dev/boiledorange73/books/caea8d4c77dbba2e23a0

次回の案内

講義予定 （全13回） • • • • • • • • • • • • • 1. 交通情報学入門 2. 地理情報システム(GIS)と時空間データベース 1 3. 地理情報システム(GIS)と時空間データベース 2 4. 交通データと計測 1: 基盤データ編 5. ゲスト講義1: 交通事業者とMaaS（藤垣洋平・小田急電鉄株式会社） 6. 交通データと計測 2: 動的データ編 7. 経路検索アルゴリズムと応用 8. 交通シミュレーションの技術と演習 1 9. ゲスト講義2: 交通ビッグデータ分析と活用の実務（太田恒平・株式会社トラフィックブレイン） 10. 交通シミュレーションの技術と演習 2 11. 高度化する交通サービス: ITS（Intelligent Transport Systems）・MaaS (Mobility as a Service) ・自動運転 12. データに基づいた交通政策の可能性 13. 交通情報学の未来（ディスカッション）

本日の課題 • 授業の感想、質問、要望などをコメントしてください • LMSで提出 〜5月31日（水） 23:59まで