PHPで大圏航路を描く

地球上の 2 地点間の最短距離を示す経路を大圏航路と呼ぶ。PHP プログラムを使って大圏航路を計算し、Google マップや Yahoo! JavaScript マップ上に描画する。あわせて、等角航路を計算してマップ上に描画し、各々の距離計算表を表示する。

プログラムを実行する

サンプル・プログラムの実行例

Google マップ表示

赤い曲線が大圏航路、黒い直線が等角航路である。

YOLP コンテンツジオコーダ API（Yahoo!JAPAN）は、世界のキーワードから緯度・経度を求める機能が弱い。たとえば「サンティアゴ」でキーワード検索すると、Google Geocoding API チリの首都を第一候補として返すが、YOLP コンテンツジオコーダ APIはドミニカ共和国の県都しか返さない。
そこで、地図上のマーカーをドラッグし、描画ボタンを押すことで航路を再描画する機能を採り入れた。これは Google マップでも利用できる。

サンプル・プログラムのダウンロード

ダウンロード

圧縮ファイルの内容
greatCircleSailing.php	サンプル・プログラム本体
pahooGeoCode.php	住所・緯度・経度に関わるクラス pahooGeoCode。使い方は「PHPで住所・ランドマークから最寄り駅を求める」などを参照。include_path が通ったディレクトリに配置すること。

準備：pahooGeoCode クラス

0026: class pahooGeoCode {
0027:     var $items;      //検索結果格納用
0028:     var $error;      //エラーフラグ
0029:     var $hits;       //検索ヒット件数
0030:     var $webapi; //直前に呼び出したWebAPI URL
0031:
0032:     //Google Cloud Platform APIキー
0033:     //https://cloud.google.com/maps-platform/
0034:     //※Google Maps APIを利用しないのなら登録不要
0035:     var $GOOGLE_API_KEY_1 = '**************************';   //HTTPリファラ用
0036:     var $GOOGLE_API_KEY_2 = '**************************';   //IP制限用
0037:
0038:     //Yahoo! JAPAN WebサービスアプリケーションID
0039:     //https://e.developer.yahoo.co.jp/register
0040:     //※Yahoo! JAPAN Webサービスを利用しないのなら登録不要
0041:     var $YAHOO_APPLICATION_ID = '*****************************';

地図サービスを利用するために、クラスファイル "pahooGeoCode.php" を使用する。組み込み関数 require_once を使って読めるディレクトリに配置する。ディレクトリは、設定ファイル php.ini に記述されているオプション設定 include_path に設定しておく。
クラスについては「PHP でクラスを使ってテキストの読みやすさを調べる」を参照されたい。

Yahoo! JAPAN Web サービスを利用するには Yahoo! JAPAN Web サービスアプリケーション IDが必要で、その入手方法は「Yahoo!JAPAN デベロッパーネットワーク - WebAPI の登録方法」を参照されたい。
また、地図として Google マップを利用するのであれば、Google Cloud Platform API キーが必要で、その入手方法は「Google Cloud Platform - WebAPI の登録方法」を参照されたい。

準備：地図サービス（WebAPI）の選択

0031: //地図サービス（WebAPI）の選択
0032: // 0：Google
0033: // 1：Yahoo!JAPAN
0034: define('GEOSERVICE', 1);

表示する地図は、Google マップと Yahoo!マップ（YOLP 地図）から選べる。
あらかじめ、定数 GEOSERVICE に値を設定すること。

Yahoo!マップ表示

なお、Yahoo!マップで大西洋を横断する航路を描くと、左図のように間違った航路が描かれる。仕様の問題と思われる。

Yahoo!マップ表示

左図のように大西洋を中央にもってくると、再描画せずとも正しい表示になる。

この問題は、Google マップでは発生しない。

サンプル・プログラムの流れ

WebAPI によって 2 地点の緯度・経度を求め、そこから航路計算を行う。

解説：大圏航路と等角航路

地球を球体とみなすとして、2 地点 P、Q の最短距離は、P、Q および球体の中心を通る平面と球面が交わる大円となる。これを大圏（たいけん）航路（Great circle route、大円コース）と呼ぶ。
なお、地球は完全な球体ではないので、厳密には大円と大圏航路には微妙なずれがある。今回は計算を簡便にする目的で、大円が大圏航路であるとしている。

視点を変えると、大圏航路は左図の赤線のように直線となる。

これに対し、極点に対して常に一定の角度で進むコースを等角航路（rhumb line）と呼ぶ。天体観測によって緯度は観測できても、正確な時刻を知るためのクロノメーターが無かった時代は、大圏航路より時間がかかるが、等角航路を使っていた。

メルカトル図法の世界地図では、大圏航路は曲線（赤線）になるが、等角航路（黒線）は直接となる。しかし、前述のように地球は球体であるから、実際には大圏航路の方が等角航路より短い。

計算式の導出は省くが、緯度（φ)、経度（δ）で表された出発点（φ₁，δ₁）、到着点（φ₂，δ₂）の大圏航路の軌跡は

$mimetex$

によって計算できる。
また、地球の半径を R とすると、2 地点間の大圏距離は

$mimetex$

によって計算できる。

解説：大圏航路の軌跡

0571: /**
0572:  * 2地点間の大圏航路軌跡を求める
0573:  * @param double $long_a, $lati_a  A地点の経度，緯度（世界測地系）
0574:  * @param double $long_b, $lati_b  B地点の経度，緯度（世界測地系）
0575:  * @param array  $points  軌跡の座標を格納
0576:  *                      [$n]['longitude'] 軌跡の経度（世界測地系）
0577:  *                      [$n]['latitude']  軌跡の緯度（世界測地系）
0578:  * @return int 座標の数
0579:  *
0580: */
0581: function greatCircleSailing($long_a, $lati_a, $long_b, $lati_b, &$points) {
0582:     $lati_a = deg2rad($lati_a);
0583:     $long_a = deg2rad($long_a);
0584:     $lati_b = deg2rad($lati_b);
0585:     $long_b = deg2rad($long_b);
0586:
0587:     $l1 = ($long_a >= 0) ? $long_a : 2 * pi() + $long_a;
0588:     $l2 = ($long_b >= 0) ? $long_b : 2 * pi() + $long_b;
0589:     $dd = $l2 - $l1;
0590:     $tt = 0.01;           //経度方向の増分
0591:     if ($dd < 0) {
0592:         $dd = abs($dd);
0593:         $tt = -$tt;
0594:     } else if ($dd > pi()) {
0595:         list($lati_a, $lati_b) = array($lati_b, $lati_a);
0596:         list($long_a, $long_b) = array($long_b, $long_a);
0597:         $dd = 2 * pi() - $dd;
0598:     }
0599:     $st = 0.0;
0600:     $cnt = 0;
0601:
0602:     //軌跡の計算
0603:     $latitude  = $lati_a;
0604:     $longitude = $long_a;
0605:     while ($st < $dd) {
0606:         if ($latitude  >= 0.5 * pi())    $latitude -= pi();
0607:         if ($longitude >= 1.0 * pi())    $longitude = $longitude - pi();
0608:         $points[$cnt]['latitude']  = rad2deg($latitude);
0609:         $points[$cnt]['longitude'] = rad2deg($longitude);
0610:         $longitude += $tt;
0611:         if ($longitude >= pi())     $longitude = $longitude - 2 * pi();
0612:         $latitude = (sin($lati_a) * sin($long_b - $longitude)) / (cos($lati_a) * sin($long_b - $long_a)) + (sin($lati_b) * sin($long_a - $longitude)) / (cos($lati_b) * sin($long_a - $long_b));
0613:         $latitude = atan($latitude);
0614:         if (sin($lati_a) / sin($lati_b) < 0)    $latitude += pi();
0615:
0616:         $st += abs($tt);
0617:         $cnt++;
0618:     }
0619:     $points[$cnt]['latitude']  = rad2deg($lati_b);
0620:     $points[$cnt]['longitude'] = rad2deg($long_b);
0621:     $cnt++;
0622:
0623:     return $cnt;
0624: }

メソッド greatCircleSailing は、前述の計算式に則り、大圏航路の軌跡を座標配列として取得する。経度方向に 0.01 ラジアンずつ増分させ、前述の計算式により緯度を計算してゆくものである。

なお、「PHP で弾道ミサイルの軌道を計算する」で紹介したように、Google マップでは、メソッド Polyline のプロパティ geodesic を true にすることで、大圏航法による最短距離を近似描画できる。しかし、Yahoo! JavaScript マップに、この機能がないため、PHP 側でメソッドを用意した。Google マップ描画時にも同じメソッドを使って軌跡を計算している。

解説：大圏航路の距離

0551: /**
0552:  * 2地点間の大圏航路距離を求める
0553:  * @param double $long_a, $lati_a  A地点の経度，緯度（世界測地系）
0554:  * @param double $long_b, $lati_b  B地点の経度，緯度（世界測地系）
0555:  * @return double 大圏航路距離
0556:  *
0557: */
0558: function greatCircleDistance($long_a, $lati_a, $long_b, $lati_b) {
0559:     $lati_a = deg2rad($lati_a);
0560:     $long_a = deg2rad($long_a);
0561:     $lati_b = deg2rad($lati_b);
0562:     $long_b = deg2rad($long_b);
0563:
0564:     //距離の計算
0565:     $ll = abs($long_b - $long_a);
0566:     $distance = 6371.0 * acos(sin($lati_a) * sin($lati_b) + cos($lati_a) * cos($lati_b) * cos($ll));
0567:
0568:     return $distance;
0569: }

メソッド greatCircleDistance は、前述の計算式に則り、大圏航路の距離を求める。

解説：等角航路の軌跡

0626: /**
0627:  * 2地点間の等角航路軌跡を求める
0628:  * @param double $long_a, $lati_a  A地点の経度，緯度（世界測地系）
0629:  * @param double $long_b, $lati_b  B地点の経度，緯度（世界測地系）
0630:  * @param array  $points  軌跡の座標を格納
0631:  *                      [$n]['longitude'] 軌跡の経度（世界測地系）
0632:  *                      [$n]['latitude']  軌跡の緯度（世界測地系）
0633:  * @return double $distance 2地点間の等角航路距離を格納
0634:  * @return int 座標の数
0635:  *
0636: */
0637: function rhumbLine($long_a, $lati_a, $long_b, $lati_b, &$points, &$distance) {
0638:     $lati_a = deg2rad($lati_a);
0639:     $long_a = deg2rad($long_a);
0640:     $lati_b = deg2rad($lati_b);
0641:     $long_b = deg2rad($long_b);
0642:     $n = 200;
0643:
0644:     //経度方向の増分
0645:     $l1 = ($long_a >= 0) ? $long_a : 2 * pi() + $long_a;
0646:     $l2 = ($long_b >= 0) ? $long_b : 2 * pi() + $long_b;
0647:     $d1 = $l2 - $l1;
0648:     $t1 = 0.01;
0649:     $n = abs($d1) / $t1;
0650:     if ($d1 < 0) {
0651:         $d1 = 2 * pi() + $d1;
0652:         $t1 = -$t1;
0653:     } else if ($d1 > pi()) {
0654:         list($lati_a, $lati_b) = array($lati_b, $lati_a);
0655:         list($long_a, $long_b) = array($long_b, $long_a);
0656:         $d1 = 2 * pi() - $d1;
0657:         $n = abs($d1) / $t1;
0658:     }
0659:
0660:     //緯度方向の増分
0661:     $l1 = $lati_a;
0662:     $l2 = $lati_b;
0663:     $d2 = $l2 - $l1;
0664:     if ($d2 >= pi()) {
0665:         $d2 -= pi();
0666:         $t2 = - $d2 / $n;
0667:     } else {
0668:         $t2 = $d2 / $n;              //緯度方向の増分
0669:     }
0670:
0671:     //軌跡の計算
0672:     $latitude  = $lati_a;
0673:     $longitude = $long_a;
0674:     $distance  = 0.0;
0675:     for ($cnt = 0; $cnt < $n; $cnt++) {
0676:         if ($latitude  >= 0.5 * pi())    $latitude -= pi();
0677:         if ($longitude >= 1.0 * pi())    $longitude = $longitude - 2 * pi();
0678:         $points[$cnt]['latitude']  = rad2deg($latitude);
0679:         $points[$cnt]['longitude'] = rad2deg($longitude);
0680:         $longitude += $t1;
0681:         $latitude  += $t2;
0682:         if ($cnt >= 1) {
0683:             $distance += $this->greatCircleDistance($points[$cnt - 1]['longitude'], $points[$cnt - 1]['latitude'], $points[$cnt]['longitude'], $points[$cnt]['latitude']);
0684:         }
0685:     }
0686:     $points[$cnt]['latitude']  = rad2deg($lati_b);
0687:     $points[$cnt]['longitude'] = rad2deg($long_b);
0688:     $distance += $this->greatCircleDistance($points[$cnt - 1]['longitude'], $points[$cnt - 1]['latitude'], $points[$cnt]['longitude'], $points[$cnt]['latitude']);
0689:     $cnt++;
0690:
0691:     return $cnt;
0692: }

メソッド rhumbLine は、等角航路の軌跡を座標配列として取得する。経度方向に 0.01 ラジアンずつ増分させ、それに応じて緯度方向もリニアに増分させてゆく。

なお、等角航路は複雑な曲線を描くため、距離を求める適当な方程式がない。そこで、隣り合う座標間の大圏距離を積分し、等角航路の距離を同時に算出する。

解説：初期値

0036: //マップの表示サイズ（単位：ピクセル）
0037: define('MAP_WIDTH',  600);
0038: define('MAP_HEIGHT', 400);
0039: //マップID
0040: define('MAPID', 'map_id');
0041: //初期値
0042: define('DEF_LONGITUDE0', 139.766667);      //地図中心（経度）
0043: define('DEF_LATITUDE0',  35.681111);      //　　　　（緯度）
0044: define('DEF_FROM',     '東京駅');            //from（クエリ）
0045: define('DEF_TO',     'ワシントン');        //to　（クエリ）
0046: define('DEF_CAT_FROM', 'landmark');        //from（カテゴリ）
0047: define('DEF_CAT_TO', 'world');            //to　（カテゴリ）
0048: define('DEF_LONGITUDE1', 139.766667);      //from（経度）
0049: define('DEF_LATITUDE1',  35.681111);      //　　（緯度）
0050: define('DEF_LONGITUDE2', -77.0329165);      //to　（経度）
0051: define('DEF_LATITUDE2',  38.8878665);     //    （緯度）
0052: define('DEF_TYPE',     'roadmap');            //マップタイプ

各種のデフォルト・パラメータは、これらの定数によって設定されている。自由に変更できる。

解説：マーカー・ドラッグ用スクリプト

0180: /**
0181:  * マーカー・ドラッグ用スクリプト：Yahoo! JavaScriptマップ
0182:  * @param bool $drag TRUE：ドラッグで再描画する／FALSE：しない（デフォルト）
0183:  * @return string JavaScript
0184: */
0185: function jsDragMarker_YOLPmap($drag=FALSE) {
0186:     $submit = $drag ? 'document.myform.submit();' : '';
0187:
0188: $js =<<< EOT
0189:     marker_A.setDraggable(true);
0190:     marker_A.bind('dragend', function() {
0191:         marker_A.setIcon(icon_A);
0192:         marker_A.setDraggable(true);
0193:         latlng = marker_A.getLatLng();
0194:         document.getElementById('longitude1').value = latlng.lng();
0195:         document.getElementById('latitude1').value  = latlng.lat();
0196:         document.getElementById('from').value = '';
0197:         {$submit}
0198:     });
0199:
0200:     marker_B.setDraggable(true);
0201:     marker_B.bind('dragend', function() {
0202:         marker_B.setIcon(icon_B);
0203:         marker_B.setDraggable(true);
0204:         latlng = marker_B.getLatLng();
0205:         document.getElementById('longitude2').value = latlng.lng();
0206:         document.getElementById('latitude2').value  = latlng.lat();
0207:         document.getElementById('to').value = '';
0208:         {$submit}
0209:     });
0210:
0211: EOT;
0212:     return $js;
0213: }

ユーザー関数 jsDragMarker_YOLPmap は、Yahoo! JavaScript マップ上のマーカーをドラッグし、別の地点に移動する JavaScript を生成する。
まず、2 つのマーカーについて、setDraggable メソッドを使ってドラッグ可能になるよう設定する。次に、ドラッグ終了時に発生する dragend イベントをフックし、マーカーを再描画するとともに、終了時の座標を取得する。
なお、引数 $drag=TRUE に設定すると、ドラッグ終了時に航路を再描画する。サーに負荷をかけないように、デフォルトでは FALSE にしている。

0244: /**
0245:  * マーカー・ドラッグ用スクリプト：Googleマップ
0246:  * @param bool $drag TRUE：ドラッグで再描画する／FALSE：しない（デフォルト）
0247:  * @return string JavaScript
0248: */
0249: function jsDragMarker_Gmap($drag=FALSE) {
0250:     $submit = $drag ? 'document.myform.submit();' : '';
0251:
0252: $js =<<< EOT
0253:     marker_A.setDraggable(true);
0254:     marker_A.addListener('dragend', function() {
0255:         marker_A.setIcon(icon_A);
0256:         marker_A.setDraggable(true);
0257:         latlng = marker_A.getPosition();
0258:         document.getElementById('longitude1').value = latlng.lng();
0259:         document.getElementById('latitude1').value  = latlng.lat();
0260:         document.getElementById('from').value = '';
0261:         {$submit}
0262:     });
0263:
0264:     marker_B.setDraggable(true);
0265:     marker_B.addListener('dragend', function() {
0266:         marker_B.setIcon(icon_B);
0267:         marker_B.setDraggable(true);
0268:         latlng = marker_B.getPosition();
0269:         document.getElementById('longitude2').value = latlng.lng();
0270:         document.getElementById('latitude2').value  = latlng.lat();
0271:         document.getElementById('to').value = '';
0272:         {$submit}
0273:     });
0274:
0275: EOT;
0276:     return $js;
0277: }

ユーザー関数 jsDragMarker_Gmap は、Google マップ上のマーカーをドラッグし、別の地点に移動する JavaScript を生成するもので、前述の jsDragMarker_YOLPmap と同じ処理。Google マップに合わせてメソッドなどを変更している。

解説：航路描画スクリプト

0215: /**
0216:  * 航路描画スクリプト：Yahoo! JavaScriptマップ
0217:  * @param int    $id     識別ID
0218:  * @param array  $items  航路の座標
0219:  * @param string $color  描画色（RGB）
0220:  * @param int    $weight 線の太さ（省略時=1）
0221:  * @return string JavaScript
0222: */
0223: function jsDrawLine_YOLPmap($id, $items, $color, $weight=1) {
0224: $js =<<< EOT
0225:     var latlng{$id} = [
0226:
0227: EOT;
0228:     foreach ($items as $item) {
0229: $js .=<<< EOT
0230:         new Y.LatLng({$item['latitude']}, {$item['longitude']}),
0231:
0232: EOT;
0233:     }
0234: $js .=<<< EOT
0235:     ];
0236:     var style{$id} = new Y.Style("{$color}", {$weight}, 1);
0237:     var polyline{$id} = new Y.Polyline(latlng{$id}, {strokeStyle: style{$id}});
0238:     ymap.addFeature(polyline{$id});
0239:
0240: EOT;
0241:     return $js;
0242: }

ユーザー関数 jsDrawLine_YOLPmap は、Yahoo! JavaScript マップ上で、計算した大圏航路または等角高度の座標を繋ぐ直線を描く JavaScript を生成する。

0279: /**
0280:  * 航路描画スクリプト：Yahoo! JavaScriptマップ
0281:  * @param int    $id     識別ID
0282:  * @param array  $items  航路の座標
0283:  * @param string $color  描画色（RGB）
0284:  * @param int    $weight 線の太さ（省略時=1）
0285:  * @return string JavaScript
0286: */
0287: function jsDrawLine_Gmap($id, $items, $color, $weight=1) {
0288: $js =<<< EOT
0289:     var pline{$id} = new google.maps.Polyline({
0290:         map: map,
0291:         path:[
0292:
0293: EOT;
0294:     foreach ($items as $item) {
0295: $js .=<<< EOT
0296:         new google.maps.LatLng({$item['latitude']}, {$item['longitude']}),
0297:
0298: EOT;
0299:     }
0300: $js .=<<< EOT
0301:     ],
0302:     strokeColor: "{$color}",
0303:     strokeOpacity: 1,
0304:     strokeWeight: {$weight},
0305:     zIndex: 1
0306:     });
0307:
0308: EOT;
0309:     return $js;
0310: }

ユーザー関数 jsDrawLine_Gmap は、Google マップ上で、計算した大圏航路または等角高度の座標を繋ぐ直線を描く JavaScript を生成するもので、前述の jsDrawLine_YOLPmap と同じ処理。

解説：距離計算表

0312: /**
0313:  * 距離計算表を作成
0314:  * @param array  $items 地点・距離情報
0315:  * @param   object $pgc   pahooGeoCodeクラス
0316:  * @return string 距離計算表（HTML）
0317: */
0318: function makeDistanceTable($items, $pgc) {
0319:     $d11 = sprintf('%.1f', round($items['dist1'], 1));
0320:     $d12 = sprintf('%.1f', round($pgc->km2mi($items['dist1']), 1));
0321:     $d13 = sprintf('%.1f', round($pgc->km2nm($items['dist1']), 1));
0322:     $d21 = sprintf('%.1f', round($items['dist2'], 1));
0323:     $d22 = sprintf('%.1f', round($pgc->km2mi($items['dist2']), 1));
0324:     $d23 = sprintf('%.1f', round($pgc->km2nm($items['dist2']), 1));
0325:     $ww  = MAP_WIDTH / 4;
0326:
0327: $html =<<< EOT
0328: <table class="plists">
0329: <caption>距離計算</caption>
0330: <tr>
0331: <th style="width:{$ww}px;"> </th>
0332: <th style="width:{$ww}px;">キロメートル</th>
0333: <th style="width:{$ww}px;">マイル</th>
0334: <th style="width:{$ww}px;">海里</th>
0335: </tr>
0336: <tr>
0337: <th>大圏航路</th>
0338: <td>{$d11}</th>
0339: <td>{$d12}</th>
0340: <td>{$d13}</th>
0341: </tr>
0342: <tr>
0343: <th>等角航路</th>
0344: <td>{$d21}</th>
0345: <td>{$d22}</th>
0346: <td>{$d23}</th>
0347: </tr>
0348: </table>
0349:
0350: EOT;
0351:     return $html;
0352: }