『事例で学ぶ経済・政策分析のためのGIS入門』のためのR入門：空間解析・可視化編【補足】
Keisuke Takano
はじめに
補足の概要
この補足編では，『事例で学ぶGIS入門』の内容には直接関係しないか，「より進んだ方法」として解説が省略されている一方，近年経済学などの社会科学分野の中でニーズが高まっている（と筆者が考えている）いくつかのトピックをピックアップしてまとめたものです．内容は以下の通りです．
- OpenStreetMapからのデータ取得
- 道路ネットワークに基づく空間解析
GISオープン教材によると，OpenStreetMapは，誰もが自由に利用できる地図（空間データ）を作成するプロジェクトであり，誰もがデータを無償でダウンロードしたり，編集を自由に行ったりできます．OSMはあくまで有志によって作られたものであり，公的データや商用データに比して，地物の網羅性は完全とは言えません．しかし，（例えば）国土数値情報のようなプラットフォームにはないような空間データも利用することができます．
OpenStreetMapからのデータ取得
パッケージの読み込み
まず，作業に必要な各種パッケージを読み込みます．この演習のメインパッケージはosmdataです．osmdataパッケージは，任意の空間範囲をカバーするOSMの各種データをRを介してダウンロードするために開発されたものです．詳細については，以下の公式WEBページを参照ください．
その他のパッケージは，既に「『事例で学ぶGIS入門』のためのR入門：空間解析・可視化編」で言及したものばかりですので，ここでの説明は省略します．
library(sf)
library(tidyverse)
library(osmdata)
library(mapview)
library(openxlsx)
mapview::mapviewOptions(fgb=FALSE)福岡市域の店舗のポイントデータを取得
OSMの地物には，キーと値という2種類のタグによって属性が付与されています．キーのサブクラスとして値があり，それらの組み合わせによって，特定のカテゴリに含まれる地物を抽出することができます．どのようなキー・値が使われているかについては，OpenStreetMap WikiのMain Featuresページを参照ください．ここでは福岡市域を対象に，キーが店舗（shop）に該当するものを抽出します．osmdataのgetbb関数を用いて地名からバウンディング・ボックスを検索し，その結果をopq関数に入力すると，その地名に対応したOSMデータにアクセスできます．その後，add_osm_feature関数でキー・値で地物を指定し空間データをダウンロードした上で，osmdata_sf関数でsf形式に変換します．
#日本・福岡県・福岡市に対応するバウンディング・ボックスを取得する
shop_fuk <- osmdata::getbb(place_name="Fukuoka Fukuoka Japan") %>%
#OSMから空間データにアクセス
osmdata::opq() %>%
#キーがshopの空間データを取得
osmdata::add_osm_feature(key="shop") %>%
#sf形式に変換
osmdata::osmdata_sf()OSMデータを取得したいバウンディング・ボックスは，getbb関数で地名から取得できる他，st_bbox関数をsfオブジェクトに適用することで取得することもできます．以下では，2015年国勢調査の小地域ポリゴンデータから福岡市域のレコードを抽出したものを用いて，OSMからデータを取得しています．
#小地域ポリゴンデータを読み込み
fuk <- sf::st_read(dsn="h27ka40.shp",
#読み込み結果を表示しない
quiet=T,
#character型変数をfactor型にしない
stringsAsFactors=F,
#読み込むshapefileが従っている文字コードを指定
options="ENCODING=CP932") %>%
#変数GST_NAME（市町村名）が「福岡市」に一致するレコードのみ残す
dplyr::filter(GST_NAME=="福岡市") %>%
#OSMで用いられる座標系（WGS84）に投影変換
sf::st_transform(crs=sf::st_crs(4326))
#OSMから空間データを取得する際，引数のバウンディング・ボックスを福岡市のポリゴンデータ（fuk）から取得する
shop_fuk <- osmdata::opq(bbox=sf::st_bbox(obj=fuk)) %>%
#KEYがshopの空間データのみを残す
osmdata::add_osm_feature(key="shop") %>%
#sf形式に変換
osmdata::osmdata_sf()取得された空間データは，データ種別（ポイントデータ・ポリゴンデータ・ラインデータ）毎に分割されたリストとして返されます．OSMの中で店舗の地物はポイントデータかポリゴンデータで返されます．ポイントデータはリスト内のosm_point要素にアクセスしてそのまま読み込めばOKです．
#list形式で返された店舗の空間データから，ポイントデータに該当する要素を取り出す．
shop_fuk_point <- shop_fuk[["osm_points"]] %>%
#施設種別に関する変数shopがNAでないレコードのみ残す
dplyr::filter(!is.na(shop)) %>%
#ID・施設名・施設種別の変数のみ残す
dplyr::select(osm_id,name,shop)ポリゴンデータは，重心点を取ってポイントデータに変換します．ポリゴンデータ・ラインデータの重心点を取りたい場合は，sfのst_centroid関数を用います．
#list形式で返された店舗の空間データから，ポリゴンデータに該当する要素を取り出す．
shop_fuk_polygon <- shop_fuk[["osm_polygons"]] %>%
#施設種別に関する変数shopがNAでないレコードのみ残す
dplyr::filter(!is.na(shop)) %>%
#ID・施設名・施設種別の変数のみ残す
dplyr::select(osm_id,name,shop) %>%
#ポリゴンの重心を取る
sf::st_centroid()ポイントデータと，（重心点に変換済みの）ポリゴンデータをrbind関数を使って縦結合します．
shop_fuk_p <- rbind(shop_fuk_point,shop_fuk_polygon)
head(shop_fuk_p)## Simple feature collection with 6 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 130.3226 ymin: 33.58202 xmax: 130.4264 ymax: 33.65198
## Geodetic CRS: WGS 84
## osm_id name shop geometry
## 312503423 312503423 マリノアシティ mall POINT (130.3226 33.59501)
## 473575921 473575921 ローソン convenience POINT (130.4237 33.58691)
## 477366418 477366418 ミニストップ convenience POINT (130.3446 33.65198)
## 485707416 485707416 BOX TOWN yes POINT (130.4168 33.61967)
## 485707433 485707433 ゆめタウン博多 supermarket POINT (130.412 33.61187)
## 495981699 495981699 サニー supermarket POINT (130.4264 33.58202)データの中身を見てみると，店舗名は変数nameに，OSMの値は変数shopに入っていることがわかります．例えば，コンビニに該当するポイントデータのみ取り出してプロットすると，以下のようになります．
shop_fuk_p_conveni <- shop_fuk_p %>%
#施設種別がコンビニに該当するレコードのみ残す
dplyr::filter(shop%in%c("convenience"))
#地図上に可視化する
mapview::mapview(x=shop_fuk_p_conveni)店舗ポイントデータに住所情報を付与
ここまでで取り出したポイントデータはsfオブジェクトですので，sfパッケージの各種関数を用いた空間データ操作が可能です．例として，先ほど読み込んだ小地域ポリゴンデータとポイントデータを空間結合することにより，小地域レベルの住所を付与してみます．まず，ポリゴン・ポイントデータを投影座標系（JGD2011・平面直角座標第2系）に変換します．
#店舗ポイントデータを投影変換
shop_fuk_p <- shop_fuk_p %>%
sf::st_transform(crs=sf::st_crs(6670))
#小地域ポリゴンデータを投影変換
fuk <- fuk %>%
sf::st_transform(crs=sf::st_crs(6670))次に，小地域ポリゴンデータで必要な変数（都道府県名KEN_NAME，市町村名GST_NAME，行政区名CSS_NAME，小地域名MOJI）のみ残した上で，属性を店舗ポイントデータに空間結合します．なお，小地域ポリゴンデータは，飛び地などが存在する場合に複数のレコードに分かれているので，残した変数を用いてグループ化した上で融合を行います．
fuk <- fuk %>%
#都道府県名・市区町村名・区名・小地域名のみ残す
dplyr::select(KEN_NAME,GST_NAME,CSS_NAME,MOJI) %>%
#飛び地を結合
dplyr::group_by(KEN_NAME,GST_NAME,CSS_NAME,MOJI) %>%
dplyr::summarise() %>%
#グループ化を解除
dplyr::ungroup()
#小地域ポリゴンの属性を店舗ポイントに空間結合
shop_fuk_p <- sf::st_join(x=shop_fuk_p,y=fuk)ポイントから座標情報を取り出しテキスト形式で書き出し
「『事例で学ぶGIS入門』のためのR入門：基礎編」での演習の都合上，ここでは抽出したポイントデータをテキスト形式に変換し，結果をデータフレーム化します．
#店舗ポイントデータの座標系を地理座標系（JGD2011）に変換
shop_fuk_p <- shop_fuk_p %>%
sf::st_transform(crs=sf::st_crs(6668))
shop_fuk_p_lonlat <- shop_fuk_p %>%
#ポイントデータのX座標（経度）・Y座標（緯度）を取得し，変数として追加する
dplyr::mutate(lon=sf::st_coordinates(.)[,"X"],
lat=sf::st_coordinates(.)[,"Y"]) %>%
#geometry属性を削除する
sf::st_set_geometry(value=NULL)上で作成した位置情報data.frameをcsv形式で書き出します．なお，OSMに登録されている地物情報は日々変化し続けるので，データの抽出を行う日によってレコード数が変動します．故に，その日に取得したデータを繰り返し使いたい場合には，出力するcsvファイルの名前に日付を入れる等して，後日上書きされないようにしておく必要があります．
write.csv(x=shop_fuk_p_lonlat,
#書き出す際のファイル名を指定
file="fuk_p_lonlat.csv",
#行名は書き出さない
row.names=FALSE,
#書き出す際の文字コードはUTF-8
fileEncoding="UTF-8")write.xlsx関数を使えば，xlsx形式でも書き出せます．データ容量等の制約に応じて，お好みの形式で書き出してください．
openxlsx::write.xlsx(x=shop_fuk_p_lonlat,
#書き出す際のファイル名を指定
file="fuk_p_lonlat.xlsx")ネットワークデータを用いない経路探索
osrmパッケージを用いることで，OSMのデータを用いて経路探索が行えます．経路探索自体は，OSRMのサーバ上でC++言語を用いて行われるので，ユーザが自身で道路等の交通ネットワークをする必要もなく，高速で経路探索を行うことが可能です．詳細は，以下の公式WEBページを参照ください．
パッケージの読み込み
先ほど読み込んだパッケージに追加して，osrmを読み込みます．以下では，osrmの各種関数の適用方法について，順を追って説明していきます．
library(osrm)移動時間行列の計算
osrmTable関数は，地点間の移動距離をまとめた行列を出力するための関数です．具体例としてここでは，福岡市中央区にある書店間の移動時間行列を計算してみます．移動時間の計算に先立ち，対象となる書店のポイントデータを抽出します．
shop_fuk_p_book_chuo <- shop_fuk_p %>%
#施設種別が書店に該当するレコードのみ残す
dplyr::filter(shop%in%c("books")) %>%
#中央区にあるレコードのみ残す
dplyr::filter(CSS_NAME%in%c("中央区"))
#地図上に可視化する
mapview::mapview(x=shop_fuk_p_book_chuo)次に，osrmTable関数で移動時間行列を計算します．引数locでは地点ID・位置座標でなるデータフレームもしくはsfオブジェクト（sfオブジェクトの場合，IDには行名が用いられます），osrm.profileでは移動手段を指定します．移動手段は徒歩（foot）の他にも自転車（bike）や自動車（car）が選べます．関数を実行すると，移動時間行列，出発地，到着地を要素として持つリストが返ってきます．
osrm::osrmTable(loc=shop_fuk_p_book_chuo,osrm.profile="foot")## $durations
## 1107257430 2497508310 3015489455 4190918838 4852895811 5875272472
## 1107257430 0.0 12.2 8.1 8.0 24.2 1.7
## 2497508310 12.2 0.0 11.2 16.0 11.4 12.5
## 3015489455 8.1 11.2 0.0 6.7 21.8 8.5
## 4190918838 8.0 16.0 6.7 0.0 26.3 7.4
## 4852895811 24.2 11.4 21.8 26.3 0.0 25.0
## 5875272472 1.7 12.5 8.5 7.4 25.0 0.0
## 6817751985 39.8 27.0 37.4 41.8 15.7 40.5
## 6957421943 43.6 30.8 41.2 45.7 19.5 43.9
## 7874937285 11.6 3.7 10.9 15.4 13.4 11.9
## 9608943708 48.8 36.0 46.4 50.9 24.7 49.1
## 6817751985 6957421943 7874937285 9608943708
## 1107257430 39.8 43.6 11.6 48.8
## 2497508310 27.0 30.8 3.7 36.0
## 3015489455 37.4 41.2 10.9 46.4
## 4190918838 41.8 45.7 15.4 50.9
## 4852895811 15.7 19.5 13.4 24.7
## 5875272472 40.5 43.9 11.9 49.1
## 6817751985 0.0 4.1 28.6 9.3
## 6957421943 4.1 0.0 32.2 7.0
## 7874937285 28.6 32.2 0.0 37.4
## 9608943708 9.3 7.0 37.4 0.0
##
## $sources
## lon lat
## 1107257430 130.4008 33.59065
## 2497508310 130.3952 33.58772
## 3015489455 130.3971 33.59283
## 4190918838 130.3981 33.59325
## 4852895811 130.3891 33.58443
## 5875272472 130.4017 33.59118
## 6817751985 130.3796 33.57883
## 6957421943 130.3782 33.57713
## 7874937285 130.3971 33.58706
## 9608943708 130.3780 33.57459
##
## $destinations
## lon lat
## 1107257430 130.4008 33.59065
## 2497508310 130.3952 33.58772
## 3015489455 130.3971 33.59283
## 4190918838 130.3981 33.59325
## 4852895811 130.3891 33.58443
## 5875272472 130.4017 33.59118
## 6817751985 130.3796 33.57883
## 6957421943 130.3782 33.57713
## 7874937285 130.3971 33.58706
## 9608943708 130.3780 33.57459osrmTable関数では，出発地と到着地を別の空間データにすることも可能です．例えば，1-2番目のレコードを出発地，3番目以降のレコードを到着地として行列を計算することが可能です．その際には，引数locではなく，出発地を表す引数srcと到着地を表す引数dstに，それぞれ出発地・到着地としたい地物の空間データを指定します．
osrm::osrmTable(src=shop_fuk_p_book_chuo[1:2,],
dst=shop_fuk_p_book_chuo[-c(1:2),],
osrm.profile="foot")## $durations
## 3015489455 4190918838 4852895811 5875272472 6817751985 6957421943
## 1107257430 8.1 8 24.2 1.7 39.8 43.6
## 2497508310 11.2 16 11.4 12.5 27.0 30.8
## 7874937285 9608943708
## 1107257430 11.6 48.8
## 2497508310 3.7 36.0
##
## $sources
## lon lat
## 1107257430 130.4008 33.59065
## 2497508310 130.3952 33.58772
##
## $destinations
## lon lat
## 3015489455 130.3971 33.59283
## 4190918838 130.3981 33.59325
## 4852895811 130.3891 33.58443
## 5875272472 130.4017 33.59118
## 6817751985 130.3796 33.57883
## 6957421943 130.3782 33.57713
## 7874937285 130.3971 33.58706
## 9608943708 130.3780 33.57459経路探索
osrmRoute関数を用いることにより，地点間の最短経路を探索し，探索結果をsfオブジェクトで得ることができます．ここでは例として，1番目のレコードの書店から2番目の書店までの最短経路を探索し，オブジェクトrouteに代入した上で，地図上にプロットしてみます．引数srcには出発地の位置座標もしくはsfオブジェクト，dstには到着地の位置座標もしくはsfオブジェクトを指定します．引数returnclassには，返り値として受け取りたいデータ形式を指定します．先ほどと同様，引数osrm.profileには「foot」を指定します．
#最短経路の探索
route <- osrm::osrmRoute(src=shop_fuk_p_book_chuo[1,],
dst=shop_fuk_p_book_chuo[2,],
returnclass="sf",
osrm.profile="foot")
#地図上に可視化する
#出発地を赤でプロット
mapview::mapview(x=shop_fuk_p_book_chuo[1,],col.regions="red") +
#到着地を青でプロット
mapview::mapview(x=shop_fuk_p_book_chuo[2,],col.regions="blue") +
#探索された経路を黒でプロット
mapview::mapview(x=route,color="black")複数地点を巡回する
ある地点から出発→いくつかの場所を経由→出発地に戻る，という移動を行う際の最短経路には，osrmTrip関数を用います．こうした経路探索を行う手順は，いわゆる巡回セールスマン問題（traveling salesman problem; TSP）と呼ばれます．ここでは例として，中央区内の書店を全て回る際の最短経路を探索します．これまでと同様，引数locでは地点ID・位置座標でなるデータフレームもしくはsfオブジェクト，引数returnclassには返り値として受け取りたいデータ形式，osrm.profileでは移動手段を指定します．
#最短経路の探索
tsp <- osrm::osrmTrip(loc=shop_fuk_p_book_chuo,
returnclass="sf",
osrm.profile="foot")
#返り値の表示
tsp## [[1]]
## [[1]]$trip
## Simple feature collection with 10 features and 4 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: 130.3766 ymin: 33.57456 xmax: 130.4018 ymax: 33.59325
## Geodetic CRS: JGD2011
## start end duration distance geometry
## 1 1107257430 7874937285 11.563333 0.8600 LINESTRING (130.4008 33.590...
## 2 7874937285 6817751985 28.623333 2.1434 LINESTRING (130.3971 33.587...
## 3 6817751985 6957421943 4.086667 0.3065 LINESTRING (130.3796 33.578...
## 4 6957421943 9608943708 7.036667 0.5283 LINESTRING (130.3782 33.577...
## 5 9608943708 4852895811 24.698333 1.8523 LINESTRING (130.378 33.5745...
## 6 4852895811 2497508310 11.436667 0.8577 LINESTRING (130.3891 33.584...
## 7 2497508310 3015489455 11.188333 0.8291 LINESTRING (130.3952 33.587...
## 8 3015489455 4190918838 6.706667 0.4983 LINESTRING (130.3971 33.592...
## 9 4190918838 5875272472 7.355000 0.5517 LINESTRING (130.3981 33.593...
## 10 5875272472 1107257430 1.743333 0.1308 LINESTRING (130.4017 33.591...
##
## [[1]]$summary
## [[1]]$summary$duration
## [1] 114.4383
##
## [[1]]$summary$distance
## [1] 8.5581経路探索の結果はリスト形式で返され，1つ目の要素には経路のラインデータ，2つ目の要素には総所要時間，3つ目の要素には総移動距離が入っています．ここでは，1つ目の要素を取り出した上で，書店の位置情報とともに地図上にプロットします．
#地図上に可視化する
#中央区内の書店をプロット
mapview::mapview(x=shop_fuk_p_book_chuo) +
#探索された経路を黒でプロット
mapview::mapview(x=tsp[[1]]$trip,color="black")到達圏の探索
osrmIsochrone関数は，とある地点から到達可能な領域（等時間線）を生成するための関数です．例として，1番目のレコードの書店について等時間線を生成してみます．引数locでは地点ID・位置座標でなるデータフレームもしくはsfオブジェクト，引数returnclassには返り値として受け取りたいデータ形式，osrm.profileでは移動手段を指定します．引数breaksは，何分までの到達圏を何分刻みで出力するかを指定する引数です．出力された等時間線を，1番目の書店の位置情報と重ねてプロットしてみましょう．
#等時間線
isochrone <- osrm::osrmIsochrone(loc=shop_fuk_p_book_chuo[1,],
#徒歩0分〜10分までの到達圏を1分区切り
breaks=seq(from=0,to=10),
returnclass="sf",
osrm.profile="foot")
#等時間線をプロット
mapview::mapview(x=isochrone) +
#探索された経路を黒でプロット
mapview::mapview(x=shop_fuk_p_book_chuo[1,],col.regions="black")ネットワークデータを用いた経路探索
osrmパッケージによる経路探索は非常に便利な反面，筆者の知る限りでは，過去のOSMのネットワークデータや，OSM以外から得られるネットワークデータに基づく経路探索は行えません．ここでは，任意のネットワークデータについて経路探索を行う方法を解説します．経路探索に使用するのは，sfnetworksパッケージです．sfnetworksパッケージは，ネットワークの扱いに比較的特化したパッケージで，最短経路探索の他，到達圏の探索や最寄施設の探索，移動距離行列の計算が可能です．詳細については，以下の公式WEBページを参照ください．
sfnetworksパッケージで扱うネットワークデータは，過去のOSMから得られる2時点の道路データです．過去のOSMデータには，ohsomeパッケージを用いてアクセスすることが可能です．ohsomeは，OSMのほぼ全ての更新履歴に関するデータを配信するサービスで，道路ネットワークに限らない，OSM上の地物の新設・消滅を追うことができます．ただし，OSMの整備状況は時点・場所によるバラツキが大きいため，分析対象地域に応じて，いつまでデータを遡ることができるかを逐一吟味する必要があります．詳細については，以下の公式WEBページを参照ください．
今回の分析対象事例は，福岡市で2018年に全通した，国道3号博多バイパス（下原～多々良中学校西交差点の3.3kmが開業）です．国道3号沿線の福岡市東部地域には，福岡空港や博多港，鉄道貨物ターミナル等の物流拠点が数多く立地しており，深刻な道路渋滞が慢性化しています．博多バイパスは，そうした国道3号の渋滞を緩和し，物流を効率化する目的で建設された高規格道路です．区間の概要については，以下の画像を参照ください．
Figure 1: 博多バイパスの開業区間（出典：国土交通省九州地方整備局福岡国道事務所）
CRAN外からのパッケージのインストール
ohsomeパッケージは，R本体や各種パッケージをダウンロードするためのWEBページであるCRAN（Comprehensive
R Archive
Network）からは現状ダウンロードできません．すなわち，install.packages関数によるインストールに対応していないので，パッケージの開発元（GitHub）から直接インストールする必要があります．
GitHubからインストールを行う際には，別途remotesパッケージが必要なので，インストールが済んでいなければそちらを先にインストールします．その上で，install_github関数を用いて，GitHubからohsomeパッケージのインストールを行います．関数内で指定されているアドレスは，ohsomeパッケージのWEBページ上にも記されています．筆者のPCには，これらパッケージはインストール済みなので，Rmdファイル実行の都合上プログラムとしては実行しません．
# install.packages("remotes")
# remotes::install_github("GIScience/ohsome-r")パッケージの読み込み
演習に用いるパッケージを読み込みます．上で言及したパッケージのほか，magrittrパッケージも読み込みます．magrittrは，パイプ演算子%>%の中でさまざまなデータ操作を行うためのパッケージで，dplyrパッケージとの相性が非常に良いです．また，tidygraphは，空間ネットワークを含む様々な形式のネットワークデータを扱う際に有用なパッケージです．詳細については，以下の公式WEBページを参照ください．
library(ohsome)
library(sfnetworks)
library(magrittr)
library(tidygraph)博多バイパス全通前後の道路ラインの取得
ohsomeパッケージを用いて，博多バイパス開通前（2017年1月1日）と開通後（2019年1月1日）の道路ラインを取得します．それに先立ち，道路ラインの取得範囲となる，福岡市博多区・東区，新宮町について，小地域ポリゴンを抽出します．
#福岡県小地域ポリゴンから福岡市博多区・東区，新宮町のレコードのみを取り出す
fuk_east <- sf::st_read(dsn="h27ka40.shp",
#読み込み結果を表示しない
quiet=T,
#character型変数をfactor型にしない
stringsAsFactors=F,
#読み込むshapefileが従っている文字コードを指定
options="ENCODING=CP932") %>%
#変数GST_NAME（市町村名）が「福岡市」に一致するレコードのみ残す
dplyr::filter(CSS_NAME=="博多区"|CSS_NAME=="東区"|CSS_NAME=="新宮町") %>%
#OSMで用いられる座標系（WGS84）に投影変換
sf::st_transform(crs=sf::st_crs(4326))続いて，取得する道路ラインの種別を指定します．種別は「クエリ」という形式でohsomeに渡すので，文字列結合を行い，リクエスト形式への変換を行います．今回は例として，市区町村道レベルまでの道路データの取得をリクエストします．
#道路ラインの分類を指定するクエリ（市区町村道ぐらいまで）
filter_form_road <- paste0("highway=",c("motorway",
"motorway_link",
"trunk",
"trunk_link",
"primary",
"primary_link",
"secondary",
"secondary_link",
"tertiary",
"tertiary_link")) %>%
paste(collapse=" or ")
#クエリの表示
filter_form_road## [1] "highway=motorway or highway=motorway_link or highway=trunk or highway=trunk_link or highway=primary or highway=primary_link or highway=secondary or highway=secondary_link or highway=tertiary or highway=tertiary_link"上で作成した道路種別のリクエストを用いて，クエリ全体を作成します．引数endpointには，データを取得するためにアクセスするサーバを指定します．ある1時点の地物情報を入手したい場合には「elements/geometry」と指定しますが，ある時点までの全ての更新履歴を取得したい場合には，「elementsFullHistory/geometry」と指定してください．boundaryには，データを取得したい空間領域を指定します．sfオブジェクトの他，バウンディング・ボックスで指定することも可能です．filterには，どの地物を取得するかを記したリクエストを指定します．clipGeometryには，空間領域で地物をクリップするかを指定します．道路のような地物は，市区町村単位などではなく〇〇号線単位で作成されているケースがあるので，取得データを小さくする目的で「true」としました．
#ohsomeに送信するクエリ全体を作成
m_road <- ohsome::ohsome_query(endpoint="elements/geometry",
#データ取得領域はfuk_eastポリゴンの領域
boundary=fuk_east,
#上で作成したクエリに対応する地物を抽出
filter=filter_form_road,
#fuk_eastポリゴンで地物を切り抜く
clipGeometry="true")作成したクエリを用いて，道路ラインをダウンロードします．関数set_timeでデータ時点を，set_properties関数でメタデータの取得を行うことを指定した上で，ohsome_post関数でデータをダウンロードします．まずは，2017年1月1日のデータをダウンロードします．
#2017年1月1日時点の道路ラインを取得
road_20170101 <- m_road %>%
#データ時点を2017年1月1日に設定するクエリ
ohsome::set_time(time="2017-01-01") %>%
#メタデータの取得を指定するクエリ
ohsome::set_properties("metadata") %>%
#クエリを元にデータをダウンロード
ohsome::ohsome_post()
#データの中身
road_20170101## Simple feature collection with 1491 features and 6 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 130.2888 ymin: 33.53619 xmax: 130.4966 ymax: 33.76579
## Geodetic CRS: WGS 84
## First 10 features:
## @changesetId @lastEdit @osmId @osmType @snapshotTimestamp
## 1 12169540 2012-07-10 06:23:11 way/67792290 way 2017-01-01
## 2 38145521 2016-03-29 14:31:38 way/41575309 way 2017-01-01
## 3 21462831 2014-04-02 18:02:12 way/43557605 way 2017-01-01
## 4 32792165 2015-07-22 02:37:23 way/41176226 way 2017-01-01
## 5 15389939 2013-03-16 22:46:56 way/41177844 way 2017-01-01
## 6 15389939 2013-03-16 22:46:56 way/175150508 way 2017-01-01
## 7 32416987 2015-07-05 01:43:29 way/186154250 way 2017-01-01
## 8 32792165 2015-07-22 02:37:23 way/201674695 way 2017-01-01
## 9 15391359 2013-03-17 04:54:30 way/210507347 way 2017-01-01
## 10 38178170 2016-03-30 19:19:48 way/210507354 way 2017-01-01
## @version geometry
## 1 4 LINESTRING (130.4552 33.721...
## 2 14 LINESTRING (130.4238 33.617...
## 3 7 LINESTRING (130.4209 33.616...
## 4 8 LINESTRING (130.4214 33.592...
## 5 6 LINESTRING (130.42 33.59207...
## 6 2 LINESTRING (130.4187 33.590...
## 7 5 LINESTRING (130.4186 33.591...
## 8 7 LINESTRING (130.4096 33.597...
## 9 3 LINESTRING (130.4187 33.589...
## 10 3 LINESTRING (130.4188 33.590...現状のデータでは，geometryを除く変数名の頭に「@」がついていますが，そのままではデータ操作がしにくいので，magrittrのset_colnames関数を用いて「@」を消去します．また，取得されたデータのgeometry形式はGEOMETRYとなっていますが，これは，通常のLINESTRING・MULTILINESTRING形式の道路ラインの他に，RELATIONという形式のレコード（道路の並び順に関する変数）が混ざっているためです．変数osmTypeが「way」のレコードが道路ラインを表す値なので，まずそれだけを抽出し，sfのst_collection_extract関数を用いて，LINESTRING・MULTILINESTRING形式のレコードのみ抽出します．最後に，sfのst_cast関数を用いてMULTILINESTRING形式のレコードをLINESTRING形式に変換し，座標系を投影座標系（JGD2011・平面直角座標第2系）に変換します．
road_20170101 <- road_20170101 %>%
#列名から「@」を削除
magrittr::set_colnames(gsub("@","",colnames(.))) %>%
#変数osmTypeが「way」のレコードのみ残す
dplyr::filter(osmType=="way") %>%
#LINESTRING・MULTILINESTRINGのみ取り出し
sf::st_collection_extract(type="LINESTRING") %>%
#MULTILINESTRINGをLINESTRINGに変換
sf::st_cast(to="LINESTRING") %>%
#座標系を変換
sf::st_transform(6670)ダウンロードされた2017年の道路ラインをmapview関数を用いてプロットしてみます．ブロック内の街路や側道は網羅されていませんが，主要道は概ね網羅されていることが分かります．道路種別のクエリで，より階層が下の道路を含めれば網羅度は上がりますが，データのサイズが飛躍的に大きくなり，場合によってはサーバからダウンロードを拒否されます．
mapview::mapview(x=road_20170101)2019年1月1日のデータについても，同じ手順でダウンロードします．先ほどは，データの取得からクレンジングに至るまでの処理を2つのコードブロックに分けましたが，今回はひとまとめにして実行します．
#2019年1月1日時点の道路ラインを取得
road_20190101 <- m_road %>%
#データ時点を2019年1月1日に設定するクエリ
ohsome::set_time(time="2019-01-01") %>%
#メタデータの取得を指定するクエリ
ohsome::set_properties("metadata") %>%
#クエリを元にデータをダウンロード
ohsome::ohsome_post() %>%
#LINESTRING・MULTILINESTRINGのみ取り出し
sf::st_collection_extract(type="LINESTRING") %>%
#列名から「@」を削除
magrittr::set_colnames(gsub("@","",colnames(.))) %>%
#変数osmTypeが「way」のレコードのみ残す
dplyr::filter(osmType=="way") %>%
#MULTILINESTRINGをLINESTRINGに変換
sf::st_cast(to="LINESTRING") %>%
#座標系を変換
sf::st_transform(6670)出発地・到着地ポイントの作成
今回の最短経路探索では，出発地として新宮中央駅，到着地として福岡空港駅を用います．そのため，OSMから福岡市博多区・東区，新宮町にある鉄道駅の位置情報を取得し（こちらは今日現在の位置情報を使います），そこからさらに，福岡空港駅・新宮中央駅のポイントを取り出します．今日現在のOSMデータのダウンロードは，「OpenStreetMapからのデータ取得」セクションで説明したものと同じ要領で行えます．
#駅データ
station_20210101 <- osmdata::opq(bbox=sf::st_bbox(obj=fuk_east)) %>%
#KEYがrailway，VALUEがstationの空間データのみを残す
osmdata::add_osm_feature(key="railway",value="station") %>%
#sf形式に変換
osmdata::osmdata_sf() %>%
#ポイントデータの要素（osm_points）のみ抽出
.$osm_points %>%
#必要な変数のみ残す
dplyr::select(osm_id,name) %>%
#小地域ポリゴン内のポイントのみに絞る
sf::st_intersection(y=sf::st_geometry(fuk_east)) %>%
#座標系を変換
sf::st_transform(6670)
head(station_20210101)## Simple feature collection with 6 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -53834.71 ymin: 67790.52 xmax: -53179.63 ymax: 70240.07
## Projected CRS: JGD2011 / Japan Plane Rectangular CS II
## osm_id name geometry
## 7712171239 7712171239 箱崎九大前 POINT (-53679.79 69200.38)
## 267598864 267598864 貝塚 POINT (-53289.05 70240.07)
## 1933859936 1933859936 箱崎 POINT (-53179.63 68684.68)
## 1250743978 1250743978 福岡貨物ターミナル POINT (-53784.23 69992.32)
## 7712171237 7712171237 馬出九大病院前 POINT (-53834.71 67790.52)
## 7712171238 7712171238 箱崎宮前 POINT (-53831.5 68374.5)#駅ポイントから福岡空港駅・新宮中央駅のポイントを取得
station_20210101_od <- station_20210101 %>%
dplyr::filter(name%in%c("福岡空港","新宮中央"))
station_20210101_od## Simple feature collection with 2 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -51214.07 ymin: 66378.78 xmax: -51053.31 ymax: 79065.62
## Projected CRS: JGD2011 / Japan Plane Rectangular CS II
## osm_id name geometry
## 7712171234 7712171234 福岡空港 POINT (-51214.07 66378.78)
## 5726043051 5726043051 新宮中央 POINT (-51053.31 79065.62)ネットワークデータへの変換とデータクリーニング
sfnetworksのas_sfnetwork関数を用いて，道路ラインをsfnetworkという形式に変換します．引数xには変換したいラインデータを，directedには，ネットワークが有向グラフならばTRUE，無向グラフならばFALSEを指定します．有向グラフのネットワークとは，例えば一方通行か否かのような情報を持つネットワークのことを差しますが，今回のデータはそうした情報を持っていないのでFALSEにします．
as_sfnetwork関数で変換されたネットワークデータの中身は以下のようになっています．ネットワーク上にあるノードがPOINT形式で，エッジ（辺・枝）がLINESTRING形式で格納されている状態です．
#全ての道路をsfnetworkに変換
road_fukuoka_2017_sfn_all <- sfnetworks::as_sfnetwork(x=road_20170101,directed=FALSE)
road_fukuoka_2017_sfn_all## # A sfnetwork with 1667 nodes and 1507 edges
## #
## # CRS: EPSG:6670
## #
## # An undirected multigraph with 243 components with spatially explicit edges
## #
## # Node Data: 1,667 × 1 (active)
## # Geometry type: POINT
## # Dimension: XY
## # Bounding box: xmin: -65610.7 ymin: 59591.75 xmax: -46676.05 ymax: 84564.04
## geometry
## <POINT [m]>
## 1 (-50490.55 80107.57)
## 2 (-50494.87 80101.28)
## 3 (-53467.24 68667.28)
## 4 (-54489.3 67047.05)
## 5 (-53733.3 68473)
## 6 (-53702.89 68008.9)
## # … with 1,661 more rows
## #
## # Edge Data: 1,507 × 9
## # Geometry type: LINESTRING
## # Dimension: XY
## # Bounding box: xmin: -65942.13 ymin: 59591.75 xmax: -46676.05 ymax: 85106.09
## from to change… lastEdit osmId osmType snapshotTimestamp
## <int> <int> <dbl> <dttm> <chr> <chr> <dttm>
## 1 1 2 1.22e7 2012-07-10 06:23:11 way/… way 2017-01-01 00:00:00
## 2 3 4 3.81e7 2016-03-29 14:31:38 way/… way 2017-01-01 00:00:00
## 3 5 6 2.15e7 2014-04-02 18:02:12 way/… way 2017-01-01 00:00:00
## # … with 1,504 more rows, and 2 more variables: version <dbl>,
## # geometry <LINESTRING [m]> sfnetworksで経路探索を行う際には，出発地から到着地に至る道路ネットワークは一貫して繋がっている必要があります．ただし，ohsomeを用いて取得された道路データでは，始点・終点がどこにも繋がっていない道路（例：敷地内の歩道）が含まれていたり，道路同士を繋ぐ枝道が欠けたりしており，そのままでは経路探索に失敗してしまいます．故に，経路探索の前処理として，作成した道路ネットワークを互いに繋がっていないサブグラフに分割し，最大のものだけを残す工程を挟む必要があります．
まず，tidygraphのto_components関数（sfnetworksには含まれないので注意）を用いて，ネットワークデータをサブグラフに分割します．分割結果はリスト形式で返ってきます．例として1要素目の中身を見ると，頂点の数が18個，辺の数が17個しかないので，恐らく最大のサブグラフではなさそうです．
#道路ネットワークを（互いに繋がっていない）サブグラフに分割
road_fukuoka_2017_sfn_all_comp <- road_fukuoka_2017_sfn_all %>%
tidygraph::to_components()
#要素数
length(road_fukuoka_2017_sfn_all_comp)## [1] 243#分割結果の1要素目（結果はlist形式なので）を表示
road_fukuoka_2017_sfn_all_comp[[1]]## # A tbl_graph: 849 nodes and 924 edges
## #
## # An undirected multigraph with 1 component
## #
## # Node Data: 849 × 1 (active)
## geometry
## <POINT [m]>
## 1 (-50490.55 80107.57)
## 2 (-50494.87 80101.28)
## 3 (-53467.24 68667.28)
## 4 (-54489.3 67047.05)
## 5 (-53704.26 65802.93)
## 6 (-53798.99 65811.55)
## # … with 843 more rows
## #
## # Edge Data: 924 × 9
## from to change… lastEdit osmId osmType snapshotTimestamp
## <int> <int> <dbl> <dttm> <chr> <chr> <dttm>
## 1 1 2 1.22e7 2012-07-10 06:23:11 way/… way 2017-01-01 00:00:00
## 2 3 4 3.81e7 2016-03-29 14:31:38 way/… way 2017-01-01 00:00:00
## 3 5 6 3.28e7 2015-07-22 02:37:23 way/… way 2017-01-01 00:00:00
## # … with 921 more rows, and 2 more variables: version <dbl>,
## # geometry <LINESTRING [m]>リスト形式で返された分割結果の各要素（サブグラフ）について，lapply関数を用いて辺の数を取得します．煩雑ではありますが，
- ネットワークデータから辺のみを取り出す
- 辺の数をカウントする
- 辺をsfオブジェクトに変換した上で，辺の数を変数として追加する
という3つの工程を，1つの関数としてまとめた上で実行します．最後に，各辺が最大のサブグラフに含まれているか否かを判定する変数largestを追加します．変数largestに基づき，辺を色分けした結果をmapview関数でプロットすると，ところどころ最大のサブグラフには含まれない辺（largestが0）があることが分かります．
#各サブグラフ（listの要素）に含まれる辺の数を取得し，ひとまとめの空間データに融合した上で変数として追加する
road_fukuoka_2017_sfn_all_comp_id <- do.call(rbind,lapply(road_fukuoka_2017_sfn_all_comp,function(x){
#[x]番目のサブグラフを取り出し
comp_edge <- x %>%
#道路ネットワークから辺（道路線）のみ取り出し
tidygraph::activate(edges) %>%
#tibble形式に変換
dplyr::as_tibble()
#サブグラフの辺の数を取得
nedges <- dim(comp_edge)[1]
comp_edge_out <- comp_edge %>%
#sf形式に変換（座標系を道路データに合わせる）
sf::st_sf(crs=sf::st_crs(road_20170101)) %>%
#辺の数を変数として追加
dplyr::mutate(nedges=nedges)
return(comp_edge_out)
})) %>%
#辺の数が最大のサブグラフのレコードには1，さもなくば0を取る変数
dplyr::mutate(largest=ifelse(nedges==max(nedges),"1","0"))
#最大のサブグラフとそれ以外を別々にプロット
mapview::mapview(x=road_fukuoka_2017_sfn_all_comp_id,zcol="largest",layer.name="Connected")最大のサブグラフに含まれる（largestが1の）辺のみを残し，MULTILINESTRING形式のレコードをLINESTRING形式に変換した上で，道路ラインをsfnetwork形式に変換します．その際，後の分析で使う為に，辺の長さをsfのst_length関数で計算し，変数として加えます．
road_fukuoka_2017_sfn <- road_fukuoka_2017_sfn_all_comp_id %>%
#最大のサブグラフに含まれるレコードのみ残す
dplyr::filter(largest=="1") %>%
#MULTILINESTRINGをLINESTRINGに変換
sf::st_cast(to="LINESTRING") %>%
#辺の長さを変数として追加
dplyr::mutate(length=sf::st_length(.)) %>%
#sfnetwork形式に変換
sfnetworks::as_sfnetwork(directed=FALSE)2017年以外のデータに対しても上記の前処理が適用できるよう，一連の工程をgen_sfn関数として定義します．関数の引数には，sf形式のネットワークデータnetobjと，座標系crsを取るようにします．関数を定義できたら，2017年（念の為）・2019年の道路ラインに対して適用します．
#上記で示した手順を関数化する
gen_sfn <- function(netobj,crs=sf::st_crs(netobj)){
net_all_comp <- netobj %>%
#全ての道路を空間ネットワークに変換
sfnetworks::as_sfnetwork(directed=FALSE) %>%
#空間ネットワークを（互いに繋がっていない）サブグラフに分割
tidygraph::to_components()
net_sub <- do.call(rbind,lapply(X=net_all_comp,FUN=function(x){
#サブグラフxを取り出し
comp_edge <- x %>%
tidygraph::activate(edges) %>%
#tibble形式に変換
dplyr::as_tibble()
#サブグラフxに含まれる辺の数を取得
nedges <- dim(comp_edge)[1]
comp_edge <- comp_edge %>%
#sf形式に変換
sf::st_sf(crs=crs) %>%
#辺の数を変数として追加
dplyr::mutate(nedges=nedges)
return(comp_edge)})) %>%
#辺の数が最大のサブグラフだけ残す
dplyr::filter(nedges==max(nedges)) %>%
#不要な変数を削除
dplyr::select(-nedges) %>%
#辺の長さを計算し，変数として追加
dplyr::mutate(length=sf::st_length(.)) %>%
#サブグラフを改めて空間ネットワークに変換
sfnetworks::as_sfnetwork(directed=FALSE)
return(net_sub)
}
#定義した関数を用いて2017年の道路をsfnetworkに変換
road_fukuoka_2017_sfn <- gen_sfn(netobj=road_20170101)
#2019年の道路をsfnetworkに変換
road_fukuoka_2019_sfn <- gen_sfn(netobj=road_20190101)最短経路の探索
sfnetworksのst_network_paths関数を用いて，最短経路探索を行います．引数xにはsfnetwork形式のネットワークデータ，fromには（例えば）sf形式の出発地のポイントデータ，toには到着地のポイントデータ，weightsにはネットワークの各辺に与えられる重みを指定します．weightsは特に指定しなければ，辺の長さが用いられます．ここでは例として，出発点を福岡空港駅，到着点を新宮中央駅とする経路探索を行います．
paths_2017 <- sfnetworks::st_network_paths(x=road_fukuoka_2017_sfn,
#出発地は1つ目の駅ポイント（福岡空港駅）
from=station_20210101_od[1,],
#到着地は2つ目の駅ポイント（新宮中央駅）
to=station_20210101_od[2,],
#ネットワークの重みはエッジの長さ
weights="length")経路探索結果の中身は以下のように，2つのリストがtibbleに格納されている形になっています．要素node_pathsには最短経路上にある頂点番号，edge_pathsには辺番号が入っています．
paths_2017## # A tibble: 1 × 2
## node_paths edge_paths
## <list> <list>
## 1 <int [48]> <int [47]>paths_2017[["node_paths"]]## [[1]]
## [1] 460 435 462 490 407 408 51 50 61 62 443 449 469 470 471 472 473 550 551
## [20] 109 108 138 122 71 72 613 509 506 504 505 521 597 308 283 293 294 297 292
## [39] 240 241 233 211 214 247 281 280 272 273paths_2017[["edge_paths"]]## [[1]]
## [1] 423 426 462 463 379 460 38 458 44 404 408 434 433 436 435 437 525 526 527
## [20] 73 93 94 84 50 619 610 481 478 477 493 585 583 274 273 258 260 271 257
## [39] 200 202 195 181 206 244 242 241 231sfnetworksで探索された最短経路をsf形式にするには一手間かかります．まず，エッジ番号リストを取り出した上で，dplyrのas_tibble関数でtibble形式に変換します．すると，data.frameに似た形式にリストが変換されます．
#エッジ番号リストを取り出し
edge_list_2017 <- paths_2017[["edge_paths"]][[1]] %>%
#tibble形式に変換
dplyr::as_tibble()
edge_list_2017## # A tibble: 47 × 1
## value
## <int>
## 1 423
## 2 426
## 3 462
## 4 463
## 5 379
## 6 460
## 7 38
## 8 458
## 9 44
## 10 404
## # … with 37 more rows
## # ℹ Use `print(n = ...)` to see more rows上で作成した辺番号リストに基づき，sfnetwork形式のネットワークデータから，最短経路に該当する辺を取り出します．まず，sfnetworksのactivate関数を用いて，sfnetwork内の要素のうち，辺の要素のみを取り出します．その後，辺要素をtibble形式に変換し（これで，辺要素がラインデータとして認識されます），上で作成した最短経路辺リストを内部結合します．最後に，探索された最短経路をmapview関数で可視化しましょう．
paths_2017 <- road_fukuoka_2017_sfn %>%
#sfnetwork内の要素のうち，エッジの要素のみを取り出し
sfnetworks::activate(edges) %>%
#tibble形式に変換
dplyr::as_tibble() %>%
#通し番号変数（エッジ番号に相当）を作成
dplyr::mutate(edge_id=1:nrow(.)) %>%
#最短経路エッジリストを内部結合
dplyr::inner_join(y=edge_list_2017,by=c("edge_id"="value"))
mapview::mapview(x=paths_2017)2017年以外の道路ネットワークに対しても，経路探索から経路のラインデータ化に至る一連の工程を適用できるよう，処理内容をst_network_paths_line関数として定義します．引数は，st_network_pathsのものと共通です．関数を定義できたら，2017年（念の為）・2019年の道路ネットワークに対して適用します．
st_network_paths_line <- function(x,from,to,weights="length"){
paths <- sfnetworks::st_network_paths(x=x,from=from,to=to,weights=weights)
#エッジ番号リストを取り出し
edge_list <- paths[["edge_paths"]][[1]] %>%
#tibble形式に変換
dplyr::as_tibble()
paths_sf <- x %>%
#sfnetwork内の要素のうち，エッジの要素のみを取り出し
sfnetworks::activate(edges) %>%
#tibble形式に変換
dplyr::as_tibble() %>%
#通し番号変数（エッジ番号に相当）を作成
dplyr::mutate(edge_id=1:nrow(.)) %>%
#最短経路エッジリストを内部結合
dplyr::inner_join(y=edge_list,by=c("edge_id"="value"))
return(paths_sf)
}
#2017年の道路ネットワークで経路探索
paths_2017 <- st_network_paths_line(x=road_fukuoka_2017_sfn,
from=station_20210101_od[1,],
to=station_20210101_od[2,],
weights="length")
#2019年の道路ネットワークで経路探索
paths_2019 <- st_network_paths_line(x=road_fukuoka_2019_sfn,
from=station_20210101_od[1,],
to=station_20210101_od[2,],
weights="length")探索結果をmapview関数を用いてプロットすると，2019年の最短経路は博多バイパス経由のものに変わっています．パイパス全通前後で最短経路の総距離を比べると，2km弱の距離短縮が実現したことが分かります．
#2017年の最短経路を青でプロット
mapview::mapview(x=paths_2017,color="blue") +
#2019年の最短経路を赤でプロット
mapview::mapview(x=paths_2019,color="red")#2017年経路の総距離
sum(paths_2017$length)## 20671.13 [m]#2019年経路の総距離
sum(paths_2019$length)## 18885.58 [m]複数の出発地・到着地ポイントについて一括経路探索
st_network_paths関数（とそれをベースに定義したst_network_paths_line関数）では，1つの出発地・到着地ポイントの組み合わせに対してのみ最短経路が探索可能です．ここでは関数を，複数の出発地・到着地ポイントについて一括で経路探索可能なように拡張します．引数from_pointsには1つもしくは複数の出発地ポイント，to_pointsには到着地ポイントを取ります．ただし，出発地・到着地ポイントに付されたID変数の名前は必ず「id」としてください．
関数が定義できたら，先ほど取得した福岡東部の駅ポイントに対して適用します．1・11個目のポイントを出発地に，21・31個目のポイントを到着地に使います．経路探索に先立ち，地物idの変数名「osm_id」を「id」に変更します．実行結果には，経路のラインデータのgeometry情報の他，出発地ポイントのidに対応するfrom_id，到着地ポイントのidに対応するto_idが含まれます．
#複数の出発地・到着地ポイントについて一括で経路探索する関数を定義
#出発地・到着地ポイントに付されたID変数の名前は必ず「id」とすること
st_network_paths_line_multi <- function(x,from_points,to_points,weights="length"){
#出発地ポイントidと到着地のポイントidの組み合わせで成るdata.frameを作成
od_combn <- expand.grid(from_points$id,to_points$id,stringsAsFactors=FALSE) %>%
#1列目の名前をfrom，2列目の名前をtoに変更
magrittr::set_colnames(c("from","to"))
#od_combnの各行について以下の関数を適用
paths_combn <- apply(X=od_combn,MARGIN=1,FUN=function(y){
#出発地ポイント（from_points）のうち，変数idの値がod_combnの1列目（from）の値と一致するものを抽出
from_point <- dplyr::filter(from_points,id==y[1])
#到着地ポイント（to_points）のうち，変数idの値がod_combnの2列目（to）の値と一致するものを抽出
to_point <- dplyr::filter(to_points,id==y[2])
#抽出された出発地・到着地ポイントに，先ほど定義したst_network_paths_line関数を適用
paths <- st_network_paths_line(x=x,from=from_point,to=to_point,weights=weights) %>%
#最短経路ラインを融合
dplyr::mutate(uni=1) %>%
dplyr::group_by(uni) %>%
dplyr::summarise() %>%
dplyr::ungroup() %>%
dplyr::select(-uni) %>%
#融合されたラインに出発地・到着地idを変数として追加
dplyr::mutate(from_id=from_point$id,
to_id=to_point$id)
return(paths)
}) %>%
#処理結果はlist形式で返ってくるので，行方向に結合
do.call(rbind,.)
return(paths_combn)
}
#出発地ポイントを作成
from_points <- station_20210101[c(1,11),] %>%
#ID変数の名前を関数を適用可能なように変更
dplyr::rename(id=osm_id)
#到着地ポイントを作成
to_points <- station_20210101[c(21,31),] %>%
dplyr::rename(id=osm_id)
#関数の実行
paths_2017_multi <- st_network_paths_line_multi(x=road_fukuoka_2017_sfn,
from_points=from_points,
to_points=to_points,
weights="length")
head(paths_2017_multi)## Simple feature collection with 4 features and 2 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: -54348.63 ymin: 66050.78 xmax: -51674.28 ymax: 77758.74
## Projected CRS: JGD2011 / Japan Plane Rectangular CS II
## # A tibble: 4 × 3
## geometry from_id to_id
## <MULTILINESTRING [m]> <chr> <chr>
## 1 ((-53090.27 71144.65, -53237.45 70994.9), (-51944.52 72869.03, … 771217… 5725…
## 2 ((-52260.28 72166.11, -52307.14 72188.57, -52513.42 72298, -525… 572443… 5725…
## 3 ((-53363.78 68897.55, -53324.87 68879.87, -53281.77 68862.78, -… 771217… 7733…
## 4 ((-52772.01 70326.89, -52658.52 70282.46, -52566.88 70250.18, -… 572443… 7733…変数from_idとto_idをキーに，得られた最短経路ラインに対して出発地・到着地ポイントのその他属性を結合します．その後，ラインデータの可視化用に，出発地・到着地ポイントの名称をひとまとめの文字列にした発着地変数を作成します．最後に，mapview関数で最短経路群をプロットします．
#関数の出力結果に，出発地・到着地の属性を結合
paths_2017_multi_info <- paths_2017_multi %>%
#出発地ポイントからgeometryを削除し，駅名変数の名前を変更した上で結合
dplyr::left_join(sf::st_drop_geometry(from_points) %>%
dplyr::rename(name_from=name),by=c("from_id"="id")) %>%
#到着地ポイントも同様
dplyr::left_join(sf::st_drop_geometry(to_points) %>%
dplyr::rename(name_to=name),by=c("to_id"="id")) %>%
#出発地名・到着地名をひとまとめの文字列にした発着地変数を作成
dplyr::mutate(name_combn=paste0(name_from,"_",name_to))
#最短経路をプロット
mapview::mapview(x=paths_2017_multi_info,
#経路は発着地変数で色分け
zcol="name_combn",
#レイヤ名を分かりやすくする
layer.name="Shortest Paths")