Introduction

　このWEBサイトに設置したリンクは，新しいタブ・ウィンドウで開くことを推奨します．バグの影響で，新しいタブ・ウィンドウでないと開けないリンクが多いためです．予めご容赦ください．

Objectives

　このRPubsは，RによるGIS分析の事例として，複数地点間の経路探索を行う方法をメモ書き程度にまとめたものです．具体的な事例として，「谷戸」という名前が入っている東京都内のバス停を自転車で全て回り，スタート地点に戻るような最短経路を探索する，ということをやってみます（谷戸地形と自転車が好きな筆者の単なる趣味です）．こうした経路探索を行う手順はいわゆる巡回セールスマン問題（traveling salesman problem; TSP）と呼ばれます．巡回セールスマン問題の最適化問題としての性質に関心がある方は，例えば数理システムさんのこちらのページを参照してください．
　TSPをGIS形式のデータに基づいて実行する為のパッケージの代表的なものが，OpenStreetMapの道路ネットワークデータを用いて様々な経路探索が行えるosrmパッケージです．以下では，GISデータの取得から探索された経路の可視化に至るまでの具体的なプロセスをお示しします．
　なお，このRPubsは，Rを用いたGIS分析に関する基礎的な知識を持っている読者向けのものです．RによるGIS分析については，別途チュートリアルサイトを用意してあるので，そちらを参照ください．

Exercise

Preparation

Load Packages

　まず，今回の実装例で用いるRパッケージを読み込みます．もし，インストールが済んでいないものがあれば，事前にインストールしておいてください．キーとなるパッケージは，先述のosrmと，データハンドリングに必須のパッケージdplyr，空間データの処理に必須のパッケージsf，インタラクティブな地図の可視化が可能なmapviewパッケージです．

library(dplyr)
library(osrm)
library(sf)
library(curl)
library(mapview)

Download data

　ここでは，演習に用いるshapefile群をまとめてダウンロードしておきましょう．使うデータは以下の2つです．

##バス停データ（東京都）
#ダウンロードしたzipファイルを格納する一時ファイルを設定
bs <- tempfile()
#zipファイルをダウンロード
curl::curl_download(url="https://nlftp.mlit.go.jp/ksj/gml/data/P11/P11-10/P11-10_13_GML.zip",destfile=bs)
#現在の作業ディレクトリにzipファイルを解凍
unzip(zipfile=bs,exdir=getwd())
#一時ファイルを削除
unlink(bs);rm(bs)
##標高データ
terrain <- tempfile()
curl::curl_download(url="https://nlftp.mlit.go.jp/ksj/gml/data/G04-d/G04-d-11/G04-d-11_5339-jgd_GML.zip",destfile=terrain)
unzip(zipfile=terrain,exdir=getwd())
unlink(terrain);rm(terrain)

Import & filter (or crop) shapefiles

　ダウンロード＆解凍したshapefileをst_read関数で読み込みます．なお，バス停データ内の属性テーブルの文字コードはCP932（Shift-JISと大体同じ）なので，optionでそのことを明示して読み込む必要があります（Macでコードを実行するときには特に重要）．
　バス停データについては，この時点でバス停名（P11_001）に「谷戸」という語句が含まれるもののみに絞ります．また，標高データについても，絞り込みを行った後のポイントデータがカバーしている空間範囲で切り取ります．かくなる空間範囲（X・Y座標の最大・最小値に関するベクトル形式で取得可能で，boundary box; bbと呼ばれる）を取得する為の関数は，sfのst_bbox関数です．また，取得した空間範囲に基づきshapefileを切り取る為に，sfのst_crop関数を用いています．

##バス停データ
#shapefileを読み込む
bs_yato <- sf::st_read(dsn=paste0(getwd(),"/P11-10_13_GML/P11-10_13-jgd-g_BusStop.shp"),crs=sf::st_crs(4612),stringsAsFactors=F,options="ENCODING=CP932",quiet=TRUE) %>%
  #「谷戸」が名前に含まれるバス停のみに絞る
  dplyr::filter(grepl("谷戸",P11_001)) %>%
  #座標系を投影座標系に変換
  sf::st_transform(crs=sf::st_crs(2451))
##標高データ
#shapefileを読み込む
terrain_yato <- sf::st_read(dsn="G04-d-11_5339-jgd_ElevationAndSlopeAngleFifthMesh.shp",crs=sf::st_crs(4612),stringsAsFactors=F,quiet=TRUE) %>%
  #座標系を投影座標系に変換
  sf::st_transform(crs=sf::st_crs(2451)) %>%
  #バス停データ（絞り込み済）がカバーする空間範囲でshapefileを切り取る
  sf::st_crop(sf::st_bbox(bs_yato)) %>%
  #メッシュ内の平均標高の変数だけ残す
  dplyr::select(G04d_002) %>%
  #標高変数を文字列型から数値型に変換
  dplyr::mutate(elev=as.numeric(G04d_002)) %>%
  #文字列型の方をドロップする
  dplyr::select(-G04d_002)

Routing with TSP

Routing using osrm

　このチュートリアルのメイン部分である，osrmを用いた経路探索を以下で行います．TSPを実行する為の関数は，osrmのosrmTrip関数です．今回指定する引数は3つで，巡回する場所に関するポイントデータloc，探索された経路が返される時の形式returnclass（sfかspを指定可能），巡回に使う交通手段osrm.profile（車，自転車，徒歩が指定可能）です．
　返り値の中身を見てみると，3つの要素を含むlistになることが分かります．要素tripにはsf形式の探索経路のラインデータ，要素summaryのdurationには経路の総所要時間[分]，要素distanceには総距離[km]が格納されています．結果を見るに，1日で回ることがあまり現実的ではなさそうですが，とりあえず続けます．

#TSPを実行
route_bs_yato <- osrm::osrmTrip(loc=bs_yato,returnclass="sf",osrm.profile="bike")
#返り値の中身を見る
head(route_bs_yato)
## [[1]]
## [[1]]$trip
## Simple feature collection with 34 features and 4 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -64369.54 ymin: -54737.97 xmax: -26170.94 ymax: -25558.71
## Projected CRS: JGD2000 / Japan Plane Rectangular CS IX
## First 10 features:
##    start end  duration distance                       geometry
## 1      1  34 40.011667   9.4658 LINESTRING (-39514.73 -4510...
## 2     34  24 15.805000   3.9257 LINESTRING (-33532.3 -51529...
## 3     24  31 11.056667   2.5955 LINESTRING (-32508.91 -5473...
## 4     31  25  8.910000   2.2208 LINESTRING (-32401.83 -5262...
## 5     25  33 20.703333   5.1654 LINESTRING (-32371.08 -5067...
## 6     33  26  5.728333   1.3518 LINESTRING (-34279.73 -4692...
## 7     26  32 10.768333   2.4876 LINESTRING (-35367.14 -4737...
## 8     32   5 20.810000   5.0225 LINESTRING (-36804.69 -4587...
## 9      5   4  3.643333   0.9028 LINESTRING (-33285.25 -4427...
## 10     4  20 12.303333   2.9890 LINESTRING (-33477.22 -4347...
## 
## [[1]]$summary
## [[1]]$summary$duration
## [1] 684.2517
## 
## [[1]]$summary$distance
## [1] 168.547

　ここで，探索された経路を可視化してみます．可視化には，パッケージmapviewのmapview関数を用います．可視化したいのは，経路の形状だけなので，sfのst_geometry関数で，ラインデータのgeometryのみを取り出して可視化します．

#返り値からラインデータを取り出し，トリップ（経路を構成する1本1本のライン）にIDを付与する
trip_bs_yato <- route_bs_yato[[1]]$trip %>%
  dplyr::mutate(id_trip=row_number())
#探索経路を地図上にプロットする
mapview::mapview(sf::st_geometry(trip_bs_yato))

Evaluation of searched route

　最後に，探索経路の特性を見ておきたいと思います．今回着目する特性は経路の高低差です．先ほど読み込んだ標高データと探索経路との共通部分を取ることで，トリップ（経路を構成する1本1本のライン）別・5次メッシュ別の平均標高が分かります．これを元に，トリップ毎に平均標高の標準偏差を計算し，高低差の指標として用います．ただし，標準偏差はトリップが長距離になるほど大きくなるので，トリップ長で基準化（？）します．
　基準化された高低差指標を用いて，探索経路を色分けし可視化します．一部異常に高低差指標が大きくなっている箇所がありますが，平地が多い多摩川北岸よりも，丘陵地が続く南岸（町田市や多摩市辺り）で指標の値が大きくなっていることが分かります．とはいえ，この指標も高低差の実態を捉えているかどうか微妙なので，今後改良の余地がありそうです．

#探索経路と標高データの共通部分を抽出
trip_bs_yato_terrain <- sf::st_intersection(trip_bs_yato,terrain_yato) %>%
  #geometryを消去し，data.frameに変換する
  sf::st_set_geometry(NULL)
#トリップID毎に，平均標高の標準偏差（高低差指標）を計算する
trip_bs_yato_terrain_sd <- trip_bs_yato_terrain %>%
  dplyr::group_by(id_trip) %>%
  dplyr::summarise(elev_sd=sd(elev)) %>%
  dplyr::select(id_trip,elev_sd)
#探索経路に高低差指標を結合
trip_bs_yato_rug <- trip_bs_yato %>%
  dplyr::left_join(trip_bs_yato_terrain_sd,by="id_trip") %>%
  #距離で高低差指標を基準化
  dplyr::mutate(elev_sd_dist=elev_sd/distance)
#探索経路を高低差指標で着色し可視化
mapview::mapview(trip_bs_yato_rug %>% dplyr::select(elev_sd_dist),layer.name="Ruggedness")