Mathematica を使った国土交通省の高速道路データから最短経路を求める

Key関数：Import, Flatten, Ordering, DeleteDuplicates, Apply, Map, Position, GeoDistance, Graph, Thread
Mathematica: 13.3.1

はじめに

※この記事は、2022年9月に実施した Webinar の内容をまとめたものです。

ある地点からある地点への最短経路は、グラフ理論における最短経路問題として知られており、道路や鉄道などの情報をグラフに帰着することで求めることができます。身近な応用例としては、鉄道などの経路案内があります。 Mathematica にもグラフ理論に関する組み込み関数が数多く用意されており、 Mathematica でもグラフを用意できれば簡単に経路を求めることが可能です。しかし、現状では実際の道路情報などをグラフ形式でダウンロードできるようなサービスはありません。

今回は国土数値情報 ( https://nlftp.mlit.go.jp/ksj/index.html ) の高速道路データからインターチェンジを「頂点」、各経路を「辺」としてグラフを構築し、指定したインターチェンジ間の最短経路を求めてみます。

グラフの構築

データの読み込み

Webinar 当時に使用した、2020年の高速道路のデータを使用します。
ダウンロードしたフォルダ「N06-20」にノートブックファイルを保存し、 NotebookDirectry として指定して進めていきます。
注：2021 年のデータを使用する場合は、「utf8」フォルダのデータを使用する必要があります。

In[]:=

SetDirectory[NotebookDirectory[]];

路線のデータを読み込む

路線のデータは「N06-20_HighwaySection.geojson」を使用します。

In[]:=

data1=Import["N06-20_HighwaySection.geojson","Data"]

Out[]=

TypeFeatureCollection,FeaturesTypeFeature,PropertiesN06_0012014,N06_0022014,N06_0039999,N06_004EA02_136002,N06_005Null,N06_006Null,N06_007伊豆縦貫自動車道,N06_0083,N06_0092,N06_0102,N06_011Null,GeometryLineGeoPosition

Number of points: 117

Lat bounds: {35.1,35.1}

Lon bounds: {139.,139.}

,

⋯1418⋯

,TypeFeature,Properties{N06_0012020,N06_0022020,N06_0039999,N06_004EA02_326003,

⋯3⋯

,N06_0085,N06_0091,N06_0104,N06_011Null},GeometryLineGeoPosition

Number of points: 46

Lat bounds: {34.6,34.6}

Lon bounds: {135.47997530999999,136.}

,NameN06-20_HighwaySection

Full expression not available

(

original memory size:

18.2 MB)

In[]:=

highway=data1[[2,2]];

路線の数を調べます。

In[]:=

Length[highway]

Out[]=

1420

1420個ある路線のデータの1つ目のデータを表示します。

In[]:=

highway[[1]]

Out[]=

TypeFeature,PropertiesN06_0012014,N06_0022014,N06_0039999,N06_004EA02_136002,N06_005Null,N06_006Null,N06_007伊豆縦貫自動車道,N06_0083,N06_0092,N06_0102,N06_011Null,GeometryLineGeoPosition

Number of points: 117

Lat bounds: {35.1,35.1}

Lon bounds: {139.,139.}



接合部のデータを読み込む

路線のデータは「N06-20_Joint.geojson」を使用します。

In[]:=

data2=Import["N06-20_Joint.geojson","Data"]

Out[]=

TypeFeatureCollection,FeaturesTypeFeature,PropertiesN06_0122015,N06_0132015,N06_0149999,N06_015EA03_513042,N06_016Null,N06_017Null,N06_018南相馬鹿島SIC,N06_0192,GeometryPoint[GeoPosition[{37.7153,140.919916480000012}]],

⋯2386⋯

,TypeFeature,PropertiesN06_0122020,N06_0132020,N06_0149999,

⋯1⋯

,N06_017Null,N06_018天美,N06_0191,GeometryPoint[GeoPosition[{34.593,135.539}]],NameN06-20_Joint

Full expression not available

(

original memory size:

3.7 MB)

In[]:=

joint=data2[[2,2]];

接合部の数を調べます。

In[]:=

Length[joint]

Out[]=

2388

2388個あるデータのうちの1つ目のデータを表示します。

In[]:=

joint[[1]]

Out[]=

TypeFeature,PropertiesN06_0122015,N06_0132015,N06_0149999,N06_015EA03_513042,N06_016Null,N06_017Null,N06_018南相馬鹿島SIC,N06_0192,GeometryPoint[GeoPosition[{37.7153,140.919916480000012}]]

位置情報を抽出

路線のデータの1つ目のデータを表示します。

In[]:=

highway[[1]]

Out[]=

Number of points: 117

Lat bounds: {35.1,35.1}

Lon bounds: {139.,139.}



1つ目のデータの位置情報を抽出します。

In[]:=

highway[[1,3,2,1]]

Out[]=

GeoPosition

Number of points: 117

Lat bounds: {35.1,35.1}

Lon bounds: {139.,139.}



全てのデータから位置情報を抽出

全ての路線のデータ及び全ての接合部のデータから位置情報を抽出します。

In[]:=

highwaypos=highway[[All,3,2,1]];

In[]:=

jointpos=joint[[All,3,2,1]];

抽出したデータを確認します。

In[]:=

highwaypos[[1;;5]]

Out[]=

GeoPosition

Number of points: 117

Lat bounds: {35.1,35.1}

Lon bounds: {139.,139.}

,GeoPosition

Number of points: 79

Lat bounds: {32.7,32.7}

Lon bounds: {132.,132.}

,GeoPosition

Number of points: 26

Lat bounds: {33.6,33.6}

Lon bounds: {134.,134.}

,GeoPosition

Number of points: 72

Lat bounds: {35.1,35.1}

Lon bounds: {132.33045505000001,132.37669138000001}

,GeoPosition

Number of points: 89

Lat bounds: {33.7,33.8}

Lon bounds: {131.,131.}



In[]:=

jointpos[[1;;5]]

Out[]=

{GeoPosition[{37.7153,140.919916480000012}],GeoPosition[{36.6555,137.00070761999999}],GeoPosition[{36.6001,136.902}],GeoPosition[{35.66,139.495}],GeoPosition[{33.5877,131.126}]}

辺の作成

辺のリストを作成する

路線上にある接合部の情報を取得する

Out[]=

フローチャートに従って接合部を求めます。

接合部を近い順にソートする

接合部の位置が路線のどこにあるかを求めます。

リストを並べ替える順番を求めます。

求めた順番で接合部をソートします。

隣接する接合部をペアにする

無向グラフの辺にする

全ての路線に対して実行

各路線上にある接合部を求めます。

表形式で確認します。

辺に変換します。

このエラーは路線779に重複する点があるため発生します。
エラーを回避するために、路線の位置情報から重複する点を削除します。

{}を削除し、リストを平坦化します。

平坦化したリストを確認します。

重み(距離)のリストを作成

路線内の各接合部の位置を求めます。

求まった位置をソートします。

位置をペアにします。

範囲指定で使える形式に変換します。

指定した範囲の位置情報を抽出します。

抽出した範囲の距離を求めます。

接合部間の距離が短いとメートル単位で返されることがあるため、単位をキロメートルに変換します。

値を抽出します。

注：QuantityMagnitude[dist] でも値を抽出可能できます。

全ての路線に対して実行

グラフの作成

辺と重み(距離)のリストが完成したので、グラフを作成します。

接合部の位置情報を抽出します。

頂点番号 -> {座標} の形式に変換します。

頂点の位置を指定してグラフを作成します。

オリジナルのグラフィックスと作成したグラフを比較します。オリジナルのグラフィックスは N06-20_HighwaySection.shp を使用します。ファイルを読み込んだ後、並べて表示します。

※分かっている問題点
　作成したグラフを確認すると、9 個の頂点が他のどの頂点とも接続されていませんでした。この頂点の接合部を確認したところ、路線上の位置とは一致せずに、点と点の間に存在していました。これらの頂点を自動で辺に変換することは難しいため、個別に手作業で修正する必要があります。ただし、最短経路を求める際の点に指定しなければ今回は支障がないため、このまま進めていきます。

最短経路を求めて結果をハイライトする

接合部の情報から名前の情報を抽出します。

青森ICと山口ICの番号を求めます。

最短経路を求めます。

結果をハイライトします。

距離を求めます。

インターの名前を入力とする関数を作成します。

Google Map の結果と比較してみます。
(https://www.google.com/intl/ja/permissions/geoguidelines/)

Google Map の結果では、新潟を通る経路が示されており、距離も異なります。 Google Map では新潟方面に向かう一般道を使っている箇所が存在することがわかりました。恐らく一番影響が大きいと思われる一般道路を考慮して、改めて最短経路を求めてみます。

一般道路の辺を追加する

象潟 - 遊佐比子 17.9km

あつみ温泉 - 朝日まほろば 40.8km

蟹沢 - 二ツ井白神 12.1km

小坂北 - 大館北 16.1km

新たな辺を辺のリストに追加します。

対応する距離を重みのリストに追加します。

グラフを作成します。

Google Map の距離からまだ差はありますが、かなり近い値になりました。このグラフにはまだ荒い部分があるため、細かい点を修正していくことで、Google Map と近い値を求めることができるようになると思われます。