北海道内15箇所にポケモンマンホールが設置されるので最適経路を計算してみた

この記事はFOSS4G Advent Calendar 2019の12月8日の記事です。 昨日の記事は、chomyさんのGDALでGCOM-Cの可視領域のデータからTrue Colorイメージを再構成してみたでした。

FOSS4G Advent Calendar、毎日楽しませて頂いております。 参加表明した時点ではMapboxについて何か書こうかなと思っていたんですが、やってみたいネタを思いついてしまったので路線変更し本記事を執筆するに至ります。ではどうぞ。

はじめに

11月末ごろに、北海道にポケモンマンホールが設置される旨、公式ウェブサイトや各種ニュースサイトで取り上げられました。

ポケモンマンホール『ポケふた』

１５箇所がそれなりに均等に全道に散らばっており、コンプリートはなかなか大変そうです。私の地元もなんと設置対象となっており、ポケモン好きとして、そしてジオ道へ進む身として、またアドベントカレンダーの良いネタが見つかったという事で、道内全１５箇所を巡る最適経路を計算してみたい！と思いやってみました。

『ポケふた』を位置情報データにする

位置情報を取得

公式サイトで地図情報と座標が取得出来ます。

QGIS上で１５箇所をプロット

マウスでカチカチ…

やったぜ。

CSVファイルに座標を出力

フィールド計算機で、属性に座標（経度・緯度）を追加します。＄xで経度(lon)を、＄yで緯度(lat)を取得可能。

そしてCSVでエクスポート。

['士別市', '44.1759660664951', '142.392400060057']
['比布町', '43.8750601338394', '142.471656275175']
['足寄町', '43.2442234401681', '143.546750238874']
['石狩市', '43.5859009826181', '141.423896468599']
['恵庭市', '42.8971470864568', '141.585830729499']
['洞爺湖町', '42.5661295066754', '140.819057243298']
['森町', '42.1088651156221', '140.572905228496']
['上ノ国町', '41.8084126507218', '140.095325019279']
['天塩町', '44.8845377760538', '141.747994037888']
['稚内市', '45.4170905393652', '141.676181866153']
['豊富町', '45.1033949651183', '141.778599640962']
['紋別市', '44.334669948041', '143.37218366193']
['遠軽町', '44.0218720548609', '143.497163771181']
['大空町', '43.9155408070508', '144.171387033476']
['新得町', '43.0826784443182', '142.832912172283']

最適経路計算

計算手法

さてデータの下準備が完了し、ここからが本番です。がしかし、最適経路の計算手法などは全く知らないのでとりあえずGoogleってみました。こういう、多数の地点の巡回ルートの最適化は一般に、巡回セールスマン問題と呼ばれています。

（答えが載ってる…最高…）

という訳でこちらの記事を参考にさせて頂きました。 組合せ最適化 - 典型問題 - 巡回セールスマン問題 - Qiita

コード

github.com

import csv
import numpy as np
from scipy.spatial import distance
from ortoolpy import tsp

asahikawa_airport = [142.4541427, 43.6708414]
shinchitose_airport = [141.6811815, 42.7876057]
shinkansen_hokuto = [140.6486832, 41.9046279]

start_point = shinkansen_hokuto
nodes = [start_point]

with open('./pokefuta_coordinates.csv') as f:
    reader = csv.reader(f)
    header = next(f)
    for row in reader:
        nodes.append([ row[2], row[1] ])

dist = distance.cdist(nodes, nodes)
print(tsp(nodes, dist))

pokefuta_coordinates.csvの中身は先ほどのcsvファイルと同じです。 経度・緯度の配列をnodesとして、それぞれの距離を要素数x要素数の行列distをつくり、tsp(traveling salesman problem)に渡すと結果が出力されます。

(12.651143975005297, [0, 10, 9, 8, 3, 5, 7, 6, 4, 14, 2, 13, 12, 11, 1])

最初の値は移動コストです。２つ目の返り値の配列が、移動コストを最小化する場合に通るべき経路を示します。常に1番目のnode（[0]のこと）がスタート地点かつゴール地点に設定されます。先ほどのcsvファイルで言うと士別市がスタート地点となります。この経路を図示すると以下になります。

コード自体は問題なさそうですね、ちゃんと最短経路っぽくなっています。 道民はみんな自家用ヘリコプターで移動するので海上ルートでも問題ありません。まあ車での移動だったとしても、森町と上ノ国町は同時に回るのが効率良いので、今回の分析には影響はないものとします（距離の乖離は大きそうだけれど、巡回の順番には影響がないものとする）。 というのが士別市発ルートですが、皆様のご旅行のプランに合わせていくつかのパターンをご覧ください（記載されている総移動距離は直線距離の総和による概算です）。

旭川空港発ルート

総移動距離：1267km

北海道第２位の都市旭川市、ではなく隣の東神楽町にあり、悪天候をものともしない堅牢さにより驚異の就航率99.7%を誇る我らが旭川空港から、道内１５箇所のマンホール巡りしてみませんか？

新千歳空港発ルート

総移動距離：1237km

ハブ空港として北海道の玄関口でありながら、飛行機に乗らずとも様々なテナントで楽しめる我らが新千歳空港から、１２００ｋｍドライブでマンホール制覇だ！

新幹線北斗駅発ルート

総移動距離：1239km

ついに来たぞ新幹線、東京から４時間くらいで着くらしいのでとりあえず来てみませんか、マンホール巡りしませんか？

終わりに

複数プランを示したものの、結局ルートは環状になったため当然ですが、順番が多少異なるだけでほとんど違いがありませんでした…。今回の記事は思いつきで始めてみましたが、思ったより早く記事に出来ました（トータル４時間くらい）。位置情報データの表示や加工は仕事や趣味で色々やってきましたが、今回の最適化のような分析はあまり経験がなく、実際今回のコードも、便利なライブラリに助けられただけでその計算手法などは全くわからないままです。わからない事だらけですが、ひとまずはこの記事を読まれた方に楽しんで頂けたなら幸いです。

参考サイト

組合せ最適化 - 典型問題 - 巡回セールスマン問題 - Qiita

pythonでのcsvファイルの読み込み - Qiita

QGISでラインをポイントに変換してｘｙ座標を取り出す - Qiita

基盤地図情報の水域や道路界等のXMLデータをひとつのベクターデータにする手順

はじめに

fgd.gsi.go.jp

国土地理院は、日本の国土に関する様々なベクターデータをオープンデータとして提供しています。 当然、オープンデータとはいえ出典の明記等の制約はありますが、それなりに自由に使えます。

ところが、提供されているファイル形式がJPGIS (GML)という、正直聴き馴染みのない、また使い馴染みのない形式です。 さらに悪い事に、データはメッシュ単位でのみ提供されています。まあこの点については、都道府県単位とかだとファイルサイズが大きすぎて止むを得ない感はあります…。

しかしながら提供されているデータは、水域、道路境界、はたまた建物ポリゴンや等高線など、非常に有用なデータです。 という訳で、マイナー形式でかつメッシュ単位に細かく区切られた基盤地図情報データを、①馴染みのある形式に変換して、②地物ごとに結合して、活用しやすくしてみましょう。

事前準備

基盤地図情報のダウンロード

「基本項目」からファイル選択へ

欲しいエリアを全て選択した上で、左下の「ダウンロードファイル確認へ」

全てにチェックを入れて「まとめてダウンロード」

ダウンロードしたファイルを展開すると、メッシュ単位で圧縮ファイルが得られます。

メッシュ単位の圧縮ファイルを展開すると、15項目分の.xmlファイルが得られます。

15の基本項目について

上記の.xmlファイルのファイル名を見ると、メッシュ番号のあとにAdmAreaなどの15種類の文字列が含まれている事がわかります。 ファイル名とデータの組み合わせは以下のとおりです。

ファイル名	データの種類	ジオメトリ
AdmArea	行政区域	Polygon
AdmBdry	行政界	Polyline
AdmPt	行政点	Point
BldA	建物エリア	Polygon
BldL	建物境界	Polyline
Cntr	等高線	Polyline
CommBdry	町字界線	Polyline
ElevPT	標高点	Point
GCP	基準点	Point
RailCL	軌道の中心線（ロープウェイなど）	Polyline
RdCompt	道路構成線	Polyline
RdEdg	道路縁	Polyline
WA	水域	Polygon
WL	水涯線	Polyline
WStrL	水部構造物線	Polyline

①馴染みのある形式に変換

基盤地図情報はマイナーな形式ゆえ、専用のビューアが用意されています。 がしかし、JPGIS形式はなんとQGISで読み込めます（3.x系のみのはず）。

なのでGeoJSONでもシェープファイルでも、QGISで煮るなり焼くなり好きに出来ますね。

②地物ごとに結合

ここから先は、水域データの作成を例として解説します。

全メッシュ分の水域データ（WA）をQGISに追加

「ベクタ」ー「データ管理ツール」ー「ベクタレイヤの結合」で結合します

注意点①

ここで、地物を拡大表示すると、どうみても一つの水域なのに地物としては分かれてしまっています。

これはひとつのメッシュのXMLデータ内でもさらに細かい（大体１キロくらいの）メッシュ単位のデータになっているからです。 これでは不便なので結合します。

接する地物同士を結合する

「ベクタ」ー「空間演算ツール」ー「ディゾルブ」で、接する地物同士を結合する

注意点②

ディゾルブしたベクタレイヤの地物を選択すると、何やら全ての地物を一括で選択してしまいます。

レイヤー情報を確認してみると、地物の数が1つになってしまっています。

これでは不便なので、バラバラの地物として扱えるよう処理します。

マルチパートをシングルパートにする

「ベクタ」ー「ジオメトリツール」ー「マルチパートをシングルパートに」

以上で完成です。 なお本手順により、当初は地物ごとに持っていた属性が失われています。 重要な属性は含まれていないはずですが、うまいこと復元する方法も記事にしたいです。

iOS向けGeoJSONビューア「Vanilla GIS for iOS」リリース。使い方など

はじめに

大体１ヶ月間くらいかけて、タイトルのとおり「Vanilla GIS for iOS」を開発しリリースしました。

Vanilla GIS for iOSとは

コンセプト

Vanilla GISという名前は、以前に作ったWebアプリと同じ名前です。そのコンセプトは「シンプルでユーザーがカスタム出来る地図」であり、今回のアプリに関しても同様です。というわけで、使用技術等は全然違いますが同じコンセプトという事で同じ名前にしています。単純に、これが良いネーミングだと思っている事と、他の良いネーミングが思いつかなかった、という理由もあります。

Vanilla GIS

• Kanahiro Iguchi
• Utilities
• Free
• 

apps.apple.com

機能

• GeoJSONファイルの読み込み
• 読み込んだデータの表示方法の変更（色、透過度等）
• 各種設定値等の保存・読み込み

使用技術

• Swift
• mapbox-gl-native
• ColorSlider

スライダーで色を選択するUIをつくれる「Color Slider」 - Labo288

使い方

アプリをインストールすると、「Vanilla GIS」フォルダがiPhone内に作成されます。

中身は「geojsons」フォルダと「maps」フォルダです。 前者には読み込みたいGeoJsonファイルを格納しておきます。 後者には、「Vanilla GIS」が保存したマップデータが格納されます。

GeoJSONファイルの読み込み

「Vanilla GIS」の地図画面の下部には各機能のボタンを表示するツールバーがあります。 左端のボタンが「ファイルを開く」ボタンです。

すると、先ほど説明した「geojsons」フォルダの中身が表示されます。 地図に追加したいGeoJsonファイルを選択すれば追加されます。

色、透過度等の変更

ポリゴンデータを読み込みました。 続けて色や透過度を調整したい場合は左から２番目のボタンを使います。

レイヤー一覧が表示されます。 右側のEditをタップするとレイヤー編集画面に繊維します。 なおこの画面では、レイヤーの前後位置、表示非表示の切替も可能です。

また、ジオメトリタイプが「Point」である場合に限り、「ヒートマップモード」のオンオフが可能です。

各種設定値等の保存・読み込み

マップデータを保存する場合、右下のフロッピーボタンを使います。

ファイル名を聞かれるのでお好きな名前にしてください。

無事保存出来たら、「maps」フォルダに「.vgs」ファイルとして保存されます。 以後、アプリ起動後に当該ファイルを選択する事でマップを読み込む事が出来ます。

今後の予定

2019/11/4

ver1.1がリリースされました。

• 不具合を修正しました
• 縮尺表示を有効にしました
• 「使い方（Usage）」を追加しました
• ベースマップの変更に対応しました（国土地理院、MIERUNE）

バグフィックスに加え、以下の機能追加を予定しています。

• 国土地理院タイル以外のベースマップへの切り替え
• GeoJsonファイル以外の読み込み（.shpファイルなど？）
• 「geojsons」フォルダ以外からのデータインポート
• 「.vgs」ファイルではなく汎用形式（GeoPackage？）での出力