エンジニアの荒木です。アプリで飲食店のマッピングなどで位置情報を扱う機会って多いですよね。
アプリで扱う位置情報とは、すなわち緯度、経度のことです。代表的なユースケースでは、飲食店の位置情報を地図（例Google Map）にプロットするものがあります。
しかし、マッピングだけならば特別な用意は不要です。Google Map API叩けばいいだけ。でも、もっと複雑な事をやろうとするとどうなるでしょうか？
たとえば、2点間の距離を測定したい、特定の地点から半径xxキロメートルのお店を探したいなどです。
国内ぐらいだと2点間の距離を三角関数で縮約してもそれほど問題にはならないかもしれません。しかし、2点間の距離が遠い、具体的には国を跨ぐような事例です。例えば国際空港間の距離を把握するためには、球面を仮定しての計算が必要になります。
球面を仮定しての計算・・・どうやるのでしょうか？
JavaScript上でクライアントで計算するならば、 Google Maps JavaScript API Distance and Area Functionsが使えます。
ただ、これデータベースに全部お任せすることもできます。データベースは、たとえばPostgreSQLならばPostGISプラグイン、MySQLならば8.0以上で扱うことで実現できます。
どのデータベースを使うかですが、PostgreSQLとMySQLはどちらでも良いです。ただ、PostgreSQL + PostGISプラグインの方が昔から使われている事もあり、ネット上での情報も多いですので、どちらかを選択するのであればPostgreSQLをオススメします。
今回は、PostgreSQLとPostGISプラグインを組み合わせて、空港間を飛行機の経路（大圏コース（後述）として近似した経路）を計算する方法を示した後、さらにその経路を平面地図上で可視化する方法をお話します。
環境
- データベース
- PostgreSQL >= 9.5 PostgreSQL 10.5 付属ドキュメント
- PostGIS >= 2.5 PostGIS 2.5.3dev マニュアル
- Docker for mac
- mdillon/postgis をつかいました。
概要
- 開発環境の構築（Dockerを利用した開発環境構築）
- 球面を仮定した2点間の距離（大圏コースの長さ）
- 球面上の2点間を結んだ線を平面地上に可視化する（大圏コースの平面上の可視化）
- おまけ経度180/-180を跨ぐ場合の処理
詳細
開発環境の構築
PostgreSQLに PostGISを用意します。
https://postgis.net/install/ こちらにインストール方法あります・・・が、これやらなくてもいいです。
いまは、Dockerを使えば手元で試すことできます。
PostGISをDocker使ってさわってみた、の手順をなぞりました。docker-compose.ymlのみ1行変更しました。
上記ブログと異なるのは - ./tmp:/tmp の箇所を追加したことだけです。
このtmpフォルダはデータのやりとりのために用意しました。
docker pull mdillon/postgis mkdir init echo 'CREATE EXTENSION postgis;' > init/init.sql
docker-compose.yml
version: '3'
services:
db:
container_name: postgis_container
image: mdillon/postgis
tty: true
environment:
POSTGRES_DB: sampledb
POSTGRES_USER: pguser
POSTGRES_PASSWORD: password
volumes:
- ./pgdata:/var/lib/postgresql/data
- ./init:/docker-entrypoint-initdb.d
- ./tmp:/tmp
やりとりするデータは、空港の緯度経度情報です。
空港の緯度経度情報は GADB PostgreSQL Create Tableを使いました。前記のリンクからダウンロードして ./tmp に保存します。
最後は、下記のコマンドでDcokerでコンテナを起動します。
docker-compose up -d
これで計算準備できました。
いやー、Docker素晴らしい。（昔は、インストールするだけで一苦労だった・・・。いい時代になりました。）
球面を仮定した2点間の距離
球面での2点間の距離の前に、球面を仮定したときの距離についてもう少し話します。
球面での2点間を結ぶ線分は大圏コースと呼ばれます。大圏とは、2点と地球の中心を通る平面を仮定して、地球の外周と交わる辺のことです。
Wikipediaによると以下の通りです。Wikipedia 大圏コース
球面と平面が交差している場合、共有する曲線は円となる。この内、球の中心を通る平面と球面がなす円を大円と呼ぶ。球面上に存在する2つの点を結ぶ弧のうち最短距離となるのは、弧が大円の一部である場合である。 地球における大円を特に大圏と呼ぶ。つまり、地球上の二点間を最短距離で結ぶには、大圏の一部であるルートを選択すればよいことになる。これが大圏コースである。
パブリック・ドメイン, リンク
実線はさまざまな大圏コース（破線は緯線）
球面を仮定した2点間の距離とは、上記の大圏コースの長さのことです。
これを PostgreSQLのPostGIS（いかPostGIS）で実行するには、 ST_Distance と geographyデータ型を使います。
例として 成田空港（NRT）と ニューヨークのジョン・F・ケネディ国際空港（JFK）間の大圏の距離を計算します。
その前に、空港のデータを読み込みます。上で行った開発環境構築において、 ./tmp 配下にコピーした空港データのsql gadb_postgresql_create_airports_table.sql を読み込みます。
docker exec -it postgis_container bash > psql -U pguser sampledb < /tmp/gadb_postgresql_create_airports_table.sql
次に読み込んだデータからNRTとJFKの緯度経度を取得します。
docker exec -it postgis_container psql -U pguser sampledb
Select lat_decimal, lon_decimal From airports WHERE iata_code = 'NRT';
lat_decimal | lon_decimal
-------------+-------------
35.765 | 140.386
(1 row)
sampledb=# Select lat_decimal, lon_decimal From airports WHERE iata_code = 'JFK';
lat_decimal | lon_decimal
-------------+-------------
40.64 | -73.779
(1 row)
- 成田空港（NRT）の緯度経度は 35.765, 140.386
- ジョン・F・ケネディ国際空港（JFK）の緯度経度は 40.64, -73.779
がわかりました。
緯度経度がわかれば、PostGISで大圏の距離計算するには ST_Distance をつかうだけです。
SELECT ST_Distance(ST_GeographyFromText('POINT(140.386 35.765)'), ST_GeographyFromText('POINT(-73.779 40.64)'));
st_distance ------------------ 10854382.9192822 (1 row)
補足すると、 ST_GeographyFromText はフォーマットを変換してPostGISで扱えるデータを所得するための関数です。所定のフォーマット Well-known text から球体を仮定したデータを取得することができます。
ここで出力される距離の単位はメートルですので、成田空港とジョン・F・ケネディ国際空港との距離は 10854 キロメートルであることがわかりました。
球面を仮定した計算でも特に意識せずに簡単にできました。
球面上の2点間を結んだ線を平面地上に可視化する
さて、何気なく、上で表示されていた平面上で図示されていた、2点間のゆがんだ線これを図示します。
こんな感じの絵です。ウェブ地図で大圏航路を表示する （Leaflet版） NRT-JFK
これ、最初どうやるのかなーと思って調べるとほとんど情報が出てこない。出てくるのは英語の情報ばかり。
理由は、推測になりますが、2点ありそうです。 * 日本国内は狭い。2点間を直線で結んでも無視できるから * 日本国内で遠い距離は縦方向（経度方向である）。経度方向は平面地図では歪みが少ない（同一経度であれば歪みはなし）から
日本語の情報が少ないのであれば、より有用かと思いますので、やり方を共有します！
ポイントは、 ST_Segmentize です。
これは、線分を分割する関数です。
具体的には以下の手順で使います。
- 2点間を結ぶ線分を作成する
- 線分をST_Segmentizeで分割する
- 分割した各線分の両端は座標をST_Transformで平面上に図示する
わかりやすく同一の緯度の例で示します。 具体的には下記の通りです。
- A地点が 45,0
- B地点が 45,90
2点間を結ぶ線分を作成します。LINESTRINGをつかいます。
SELECT 'SRID=4326;LINESTRING(0 45, 90 45)'::geography;
線分をST_Segmentizeで分割します。 ST_Segmentizeは対象の線分を指定した距離以下の線分で分割する関数です。
対象線分の長さは、上述の通り以下で求めることができ余す。
SELECT ST_Distance(ST_GeographyFromText('POINT(0 45)'), ST_GeographyFromText('POINT(90 45)'));
st_distance ------------------ 6690232.93254272 (1 row)
2等分したい場合は、2等分の距離を指定します。 6690232.933なので3等分すると 6690232.933 / 2 = 3345116.466 です。 なお、緯度経度を出力するために、ST_Segmentize後にST_AsTextを適用します。
SELECT ST_AsText( ST_Segmentize('LINESTRING(0 45, 90 45)'::geography, 3345116.466) );
ST_Segmentizeは位置情報と、分割する最長距離を単位をメートルで指定します。
st_astext -------------------------------------------- LINESTRING(0 45,45 54.7356103172453,90 45) (1 row)
これを平面上にプロットすると以下の通りです。
着目は中央の点の緯度が54.7356103172453であることです。
同一緯度を結んだ線分なのですが、球面を仮定した最短ルートを図辞するとこのように北半球であれば北方向に膨らみます。 これが、最初に平面で図示した時にでる歪みです。
あとは、目を細かくします。
目を細かくするには、分割する最大距離を短くするだけです。
例として20等分で描きます。
6690232.93 / 20 = 334511.647
SELECT ST_AsText( ST_Segmentize('LINESTRING(0 45, 90 45)'::geography, 334511.6466) );
これを出力すると
LINESTRING(0 45,2.20646928941875 46.0615561924978,4.49760655388603 47.0790146925829,6.87581493431254 48.0490340277217,9.34287017526721 48.9681247572473,11.8997601886993 49.8326740369412,14.546514506294 50.6389793389253 ,17.2820309129417 51.3832922176961,20.1039093610171 52.0618725728145,23.0083059648552 52.6710532092982,25.9898219857557 53.2073136423205,29.0414436490398 53.6673610767564,32.1545477812505 54.0482154060969,35.31898515120 15 54.3472940593392,38.5232478987244 54.562491753198,41.7547198768056 54.6922498596416,45 54.7356103172453,48.2452801231944 54.6922498596416,51.4767521012756 54.562491753198,54.6810148487985 54.3472940593392,57.84545221 87495 54.0482154060969,60.9585563509602 53.6673610767564,64.0101780142443 53.2073136423205,66.9916940351448 52.6710532092982,69.8960906389829 52.0618725728145,72.7179690870583 51.3832922176961,75.453485493706 50.6389793 389253,78.1002398113007 49.8326740369412,80.6571298247328 48.9681247572473,83.1241850656875 48.0490340277217,85.502393446114 47.0790146925829,87.7935307105813 46.0615561924978,90 45) (1 row)
上記の緯度経度をプロットするとExcelにて描画すると以下の通りです。
これで平面上に大圏コースを描画できました。
おまけ経度180/-180を跨ぐ場合の処理
NRT-JFK間も描画します。問題は、経度180/-180を跨ぐことです。
まずは、前述の例と同様におこないます。
SELECT ST_AsText( ST_Segmentize('LINESTRING(140.386 35.765, -73.779 40.64)'::geography, 500000) );
LINESTRING(140.386 35.765,142.056761024068 38.5020118941837,143.858958215866 41.2134179400642,145.816184688416 43.8 943065099063,147.9567216824 46.5386545342239,150.314607581444 49.1390092168574,152.930885244867 51.6860707355252,155 .854972048755 54.1681473091761,159.145984158024 56.5704518360992,162.873623704779 58.8742160203223,167.117831618251 61.0556265380494,171.965745644171 63.0846629387923,177.503611440781 64.9240771661033,-176.199405476986 66.5290385208 116,-169.11775005771 67.8483516330449,-161.303878267811 68.8284074767936,-152.918001134129 69.4205780085609,-144.230 620537686 69.5910123694109,-135.580232139591 69.3292595740188,-127.296611509064 68.6512350817493,-119.628773404718 6 7.5947605342302,-112.71253956721 66.2102974534372,-106.580860585696 64.5514438835143,-101.196202241477 62.6682644687 681,-96.4839136127109 60.604011406961,-92.3565269469807 58.3943347182076,-88.7279703058496 56.0678409994953,-85.5203 428677163 53.6471761431289,-82.6662199408836 51.1501795498224,-80.1086425638423 48.5909133088849,-77.8000880260482 4 5.9805082328162,-75.7011128926647 43.3278311192678,-73.779 40.64) (1 row)
ん？描画が変になります。
これを解決するには、JFK側の負の経度の座標を360度加えればば良いです。
この経度修正を自動的に計算してくれる関数がPostGISには用意されています。
ジオメトリの全ての地物の全てのポイント/頂点を読み、経度値が0より小さい場合には360を足します。結果は、180度を中心にした地図に描画される、経度が0-360の区間にあるデータです。
これを適用するだけです。
SELECT ST_AsText( ST_ShiftLongitude( ST_Segmentize('LINESTRING(140.386 35.765, -73.779 40.64)'::geography, 500000)::geometry ) ) ;
LINESTRING(140.386 35.765,142.056761024068 38.5020118941837,143.858958215866 41.2134179400642,145.816184688416 43.8943065099063,147.9567216824 46.5386545342239,150.314607581444 49.1390 092168574,152.930885244867 51.6860707355252,155.854972048755 54.1681473091761,159.145984158024 56.5704518360992,162.873623704779 58.8742160203223,167.117831618251 61.0556265380494,171.9 65745644171 63.0846629387923,177.503611440781 64.9240771661033,183.800594523014 66.5290385208116,190.88224994229 67.8483516330449,198.696121732189 68.8284074767936,207.081998865871 69.4 205780085609,215.769379462314 69.5910123694109,224.419767860409 69.3292595740188,232.703388490936 68.6512350817493,240.371226595282 67.5947605342302,247.28746043279 66.2102974534372,253 .419139414304 64.5514438835143,258.803797758523 62.6682644687681,263.516086387289 60.604011406961,267.643473053019 58.3943347182076,271.27202969415 56.0678409994953,274.479657132284 53. 6471761431289,277.333780059116 51.1501795498224,279.891357436158 48.5909133088849,282.199911973952 45.9805082328162,284.298887107335 43.3278311192678,286.221 40.64) (1 row)
ウェブ地図で大圏航路を表示する （Leaflet版） NRT-JFKと同じ絵が描けました！
最後に
データベースを利用した位置情報の利用例をお話しました。
位置情報久しぶりに使うと、学びがあってとても面白いですね。特に、PostGISには自分が知らない関数が山ほどあることをしれて新たな発見がありました。
調べる過程でPostGISのリファレンス読んでいると、位置情報使って作りたいものがあれこれ浮かびます。
温故知新。これ、あの案件で使ったら便利そう・・・と想像が膨らみます。
HerokuでもPostGIS使えるので、導入の敷居はありません。
位置情報を使ったサービスを作りたい！というかたは、自前でサーバーたてるという選択肢も検討してみてください！
宣伝です！
iOSとAndroidで使える、気づきの習慣化をサポートするアプリ Stockrをリリースしました！
日常のきづき、秒で忘れちゃう人 🙋 のためのアプリです。
なにも考えずアイディアをとりあえずためておくだけ。
後日その気づきを強制的に（これ大事）振り返らせてくれます。
そんな、ものぐささん 🙋 にもおすすめです。
