読者です 読者をやめる 読者になる 読者になる

y_uti のブログ

統計、機械学習、自然言語処理などに興味を持つエンジニアの技術ブログです

geospatial 拡張モジュールで二点間の距離を求める

PHP の geospatial 拡張モジュール*1を利用して、緯度と経度で表された二地点間の距離を計算してみます。
github.com

まず、拡張モジュールを取得してビルドします。通常の手順どおり phpize を利用してビルドできます*2。私の手元では、PHP 5.3.29, PHP 5.6.29, PHP 7.1.0 の各バージョンで問題なくビルドできました。

$ git clone https://github.com/php-geospatial/geospatial.git
$ cd geospatial
$ phpize && ./configure && make && make install

例として、東京からロンドンまでの距離を計算してみます。それぞれの緯度と経度は、Wikipedia によると以下のとおりでした。

Coordinates: 35°41′N 139°41′E

Tokyo - Wikipedia

Coordinates: 51°30′26″N 0°7′39″W

London - Wikipedia

次のようなコードで二点間の距離を簡単に計算できます。$tokyo, $london の定義のとおり、経度、緯度の順番で座標を指定します。データが度分秒の形式で与えられている場合は、コード例のように dms_to_decimal 関数で変換できます。vincenty 関数で二点間の距離を計算して、結果を出力します。

<?php
$tokyo = [
    'type' => 'Point',
    'coordinates' => [
        dms_to_decimal(139, 41,  0, 'E'),
        dms_to_decimal( 35, 41,  0, 'N'),
    ],
];
$london = [
    'type' => 'Point',
    'coordinates' => [
        dms_to_decimal(  0,  7, 39, 'W'),
        dms_to_decimal( 51, 30, 26, 'N'),
    ],
];
$distance = vincenty($tokyo, $london);
echo "東京からロンドンまでの距離 = " . ($distance / 1000) . "km\n";

実行結果は次のとおりです。得られた数値を国土地理院測量計算サイト での計算結果と比較したところでは、完全に一致していました*3

$ php sample.php
東京からロンドンまでの距離 = 9582.636857km

距離の計算を行う vincenty 関数は、地球を回転楕円体と見なして精度の高い計算を実行する方法です。これより精度が低いものの高速な計算法として haversine 関数も提供されています。それぞれ Wikipedia に説明されています*4。東京からロンドンまでの距離を haversine 関数で計算したところ、9569.7513961155 km と求まりました。
Vincenty法 - Wikipedia
Haversine formula - Wikipedia

geospatial モジュールには、この他にもいくつかの機能が関数として提供されています。たとえば、fraction_along_gc_line は二点間の内分点を計算する関数です。次のコード例では、東京とロンドンの間を 5 等分する座標を求めて表示します*5

<?php
for ($i = 0; $i <= 5; ++$i) {
    $point = fraction_along_gc_line($tokyo, $london, $i / 5);
    printf("北緯 %8.5f, 東経 %9.5f\n", $point['coordinates'][1], $point['coordinates'][0]);
}

実行結果は次のとおりです。最初の行は東京の座標、最後の行はロンドンの座標で、これを 5 等分するような点が出力されています。

$ php fraction_sample.php
北緯 35.68333, 東経 139.68333
北緯 50.98803, 東経 128.78428
北緯 64.32668, 東経 108.02193
北緯 70.90759, 東経  65.00509
北緯 64.72214, 東経  21.20312
北緯 51.50722, 東経  -0.12750

コード例では 5 等分した結果を示しましたが、さらに細かく内分点を計算して地図上に重ねてみると、飛行機のフライトで見るような線分を描くこともできます。下図は 500 等分した結果をプロットした例です*6
f:id:y_uti:20161225120315p:plain

なお、今回の記事では geospatial 拡張モジュールを利用しましたが、同様の機能を PHP で実装したライブラリとして phpgeo というのもあるようです。これは composer でインストールして利用できます。まだ十分に試せていないのですが、ウェブサイトの説明をざっと読んだ限りでは、サポートしている機能はこちらの方が豊富かもしれません。

*1:まだ version 0.1 という状態なので、利用する際には注意してください。

*2:PECL にも登録されているので (PECL :: Package :: geospatial)、そちらからインストールする方法もあります。

*3:リンク先のページから「No. 2 距離と方位角の計算」を辿って確認できます。

*4:Haversine 法については日本語のページは作成されていないようです。大円距離 - Wikipedia の「コンピュータによる計算」の項に書かれている数式が対応します。

*5:$tokyo, $london の定義は省略しています。

*6:計算結果を一度 csv に出力して、D3.js で描画しました。