目次
サンプル・プログラムの実行例
YOLPコンテンツジオコーダAPI（Yahoo!JAPAN）は、世界のキーワードから緯度・経度を求める機能が弱い。たとえば「サンティアゴ」でキーワード検索すると、Google Geocoding API チリの首都を第一候補として返すが、YOLPコンテンツジオコーダAPIはドミニカ共和国の県都しか返さない。
そこで、地図上のマーカーをドラッグし、描画ボタンを押すことで航路を再描画する機能を採り入れた。これは Googleマップ でも利用できる。
サンプル・プログラムのダウンロード
| greatCircleSailing.php | サンプル・プログラム本体 |
| pahooGeoCode.php | 住所・緯度・経度に関わるクラス pahooGeoCode。 使い方は「PHPで住所・ランドマークから最寄り駅を求める」などを参照。include_path が通ったディレクトリに配置すること。 |
準備：pahooGeoCode クラス
pahooGeoCode.php
37: class pahooGeoCode {
38: var $items; //検索結果格納用
39: var $error; //エラー・フラグ
40: var $errmsg; //エラー・メッセージ
41: var $hits; //検索ヒット件数
42: var $webapi; //直前に呼び出したWebAPI URL
43:
44: //Google Cloud Platform APIキー
45: //https://cloud.google.com/maps-platform/
46: //※Google Maps APIを利用しないのなら登録不要
47: var $GOOGLE_API_KEY_1 = '**************************'; //HTTPリファラ用
48: var $GOOGLE_API_KEY_2 = '**************************'; //IP制限用
49:
50: //Yahoo! JAPAN Webサービス アプリケーションID
51: //https://e.developer.yahoo.co.jp/register
52: //※Yahoo! JAPAN Webサービスを利用しないのなら登録不要
53: var $YAHOO_APPLICATION_ID = '*****************************';
クラスについては「PHPでクラスを使ってテキストの読みやすさを調べる」を参照されたい。
Yahoo! JAPAN Webサービスを利用するには Yahoo! JAPAN Webサービス アプリケーションIDが必要で、その入手方法は「Yahoo!JAPAN デベロッパーネットワーク - WebAPIの登録方法」を参照されたい。
また、地図としてGoogleマップを利用するのであれば、Google Cloud Platform APIキー が必要で、その入手方法は「Google Cloud Platform - WebAPIの登録方法」を参照されたい。
準備：地図サービスの選択
住所検索サービスは、Google、Yahoo!JAPAN、HeartRails Geo APIから選択可能で、定数 GEOSERVICE に値を設定すること。
サンプル・プログラムの流れ
解説：大圏航路と等角航路
なお、地球は完全な球体ではないので、厳密には大円と大圏航路には微妙なずれがある。今回は計算を簡便にする目的で、大円が大圏航路であるとしている。
これに対し、極点に対して常に一定の角度で進むコースを等角航路（rhumb line）と呼ぶ。天体観測によって緯度は観測できても、正確な時刻を知るためのクロノメーターが無かった時代は、大圏航路より時間がかかるが、等角航路を使っていた。
\[ \displaystyle tan(\delta) = \frac{sin(\delta_1) sin(\phi_2 - \phi)}{cos(\delta_1) sin(\phi_2 - \phi_1)} + \frac{sin(\delta_2) sin(\phi_1 - \phi)}{cos(\delta_2) sin(\phi_1 - \phi_2)} \]
によって計算できる。
また、地球の半径をRとすると、2地点間の大圏距離は
\[ \displaystyle R cos^{-1}(sin(\phi_1) sin(\delta_2) + cos(\phi_1) cos(\phi_2) cos(abs(\delta_2 - \delta_1)) \]
によって計算できる。
解説：大圏航路の軌跡
pahooGeoCode.php
712: /**
713: * 2地点間の大圏航路軌跡を求める
714: * @param float $long_a, $lati_a A地点の経度，緯度（世界測地系）
715: * @param float $long_b, $lati_b B地点の経度，緯度（世界測地系）
716: * @param array $points 軌跡の座標を格納
717: * [$n]['longitude'] 軌跡の経度（世界測地系）
718: * [$n]['latitude'] 軌跡の緯度（世界測地系）
719: * @return int 座標の数
720: *
721: */
722: function greatCircleSailing($long_a, $lati_a, $long_b, $lati_b, &$points) {
723: $lati_a = deg2rad($lati_a);
724: $long_a = deg2rad($long_a);
725: $lati_b = deg2rad($lati_b);
726: $long_b = deg2rad($long_b);
727:
728: $l1 = ($long_a >= 0) ? $long_a : 2 * pi() + $long_a;
729: $l2 = ($long_b >= 0) ? $long_b : 2 * pi() + $long_b;
730: $dd = $l2 - $l1;
731: $tt = 0.01; //経度方向の増分
732: if ($dd < 0) {
733: $dd = abs($dd);
734: $tt = -$tt;
735: } else if ($dd > pi()) {
736: list($lati_a, $lati_b) = array($lati_b, $lati_a);
737: list($long_a, $long_b) = array($long_b, $long_a);
738: $dd = 2 * pi() - $dd;
739: }
740: $st = 0.0;
741: $cnt = 0;
742:
743: //軌跡の計算
744: $latitude = $lati_a;
745: $longitude = $long_a;
746: while ($st < $dd) {
747: if ($latitude >= 0.5 * pi()) $latitude -= pi();
748: if ($longitude >= 1.0 * pi()) $longitude = $longitude - pi();
749: $points[$cnt]['latitude'] = rad2deg($latitude);
750: $points[$cnt]['longitude'] = rad2deg($longitude);
751: $longitude += $tt;
752: if ($longitude >= pi()) $longitude = $longitude - 2 * pi();
753: $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));
754: $latitude = atan($latitude);
755: if (sin($lati_a) / sin($lati_b) < 0) $latitude += pi();
756:
757: $st += abs($tt);
758: $cnt++;
759: }
760: $points[$cnt]['latitude'] = rad2deg($lati_b);
761: $points[$cnt]['longitude'] = rad2deg($long_b);
762: $cnt++;
763:
764: return $cnt;
765: }
なお、「PHPで弾道ミサイルの軌道を計算する」で紹介したように、Googleマップでは、メソッド Polyline のプロパティ geodesic を true にすることで、大圏航法による最短距離を近似描画できる。
解説：大圏航路の距離
pahooGeoCode.php
692: /**
693: * 2地点間の大圏航路距離を求める
694: * @param float $long_a, $lati_a A地点の経度，緯度（世界測地系）
695: * @param float $long_b, $lati_b B地点の経度，緯度（世界測地系）
696: * @return float 大圏航路距離
697: *
698: */
699: function greatCircleDistance($long_a, $lati_a, $long_b, $lati_b) {
700: $lati_a = deg2rad($lati_a);
701: $long_a = deg2rad($long_a);
702: $lati_b = deg2rad($lati_b);
703: $long_b = deg2rad($long_b);
704:
705: //距離の計算
706: $ll = abs($long_b - $long_a);
707: $distance = 6371.0 * acos(sin($lati_a) * sin($lati_b) + cos($lati_a) * cos($lati_b) * cos($ll));
708:
709: return $distance;
710: }
解説：等角航路の軌跡
pahooGeoCode.php
767: /**
768: * 2地点間の等角航路軌跡を求める
769: * @param float $long_a, $lati_a A地点の経度，緯度（世界測地系）
770: * @param float $long_b, $lati_b B地点の経度，緯度（世界測地系）
771: * @param array $points 軌跡の座標を格納
772: * [$n]['longitude'] 軌跡の経度（世界測地系）
773: * [$n]['latitude'] 軌跡の緯度（世界測地系）
774: * @return float $distance 2地点間の等角航路距離を格納
775: * @return int 座標の数
776: *
777: */
778: function rhumbLine($long_a, $lati_a, $long_b, $lati_b, &$points, &$distance) {
779: $lati_a = deg2rad($lati_a);
780: $long_a = deg2rad($long_a);
781: $lati_b = deg2rad($lati_b);
782: $long_b = deg2rad($long_b);
783: $n = 200;
784:
785: //経度方向の増分
786: $l1 = ($long_a >= 0) ? $long_a : 2 * pi() + $long_a;
787: $l2 = ($long_b >= 0) ? $long_b : 2 * pi() + $long_b;
788: $d1 = $l2 - $l1;
789: $t1 = 0.01;
790: $n = abs($d1) / $t1;
791: if ($d1 < 0) {
792: $d1 = 2 * pi() + $d1;
793: $t1 = -$t1;
794: } else if ($d1 > pi()) {
795: list($lati_a, $lati_b) = array($lati_b, $lati_a);
796: list($long_a, $long_b) = array($long_b, $long_a);
797: $d1 = 2 * pi() - $d1;
798: $n = abs($d1) / $t1;
799: }
800:
801: //緯度方向の増分
802: $l1 = $lati_a;
803: $l2 = $lati_b;
804: $d2 = $l2 - $l1;
805: if ($d2 >= pi()) {
806: $d2 -= pi();
807: $t2 = - $d2 / $n;
808: } else {
809: $t2 = $d2 / $n; //緯度方向の増分
810: }
811:
812: //軌跡の計算
813: $latitude = $lati_a;
814: $longitude = $long_a;
815: $distance = 0.0;
816: for ($cnt = 0; $cnt < $n; $cnt++) {
817: if ($latitude >= 0.5 * pi()) $latitude -= pi();
818: if ($longitude >= 1.0 * pi()) $longitude = $longitude - 2 * pi();
819: $points[$cnt]['latitude'] = rad2deg($latitude);
820: $points[$cnt]['longitude'] = rad2deg($longitude);
821: $longitude += $t1;
822: $latitude += $t2;
823: if ($cnt >= 1) {
824: $distance += $this->greatCircleDistance($points[$cnt - 1]['longitude'], $points[$cnt - 1]['latitude'], $points[$cnt]['longitude'], $points[$cnt]['latitude']);
825: }
826: }
827: $points[$cnt]['latitude'] = rad2deg($lati_b);
828: $points[$cnt]['longitude'] = rad2deg($long_b);
829: $distance += $this->greatCircleDistance($points[$cnt - 1]['longitude'], $points[$cnt - 1]['latitude'], $points[$cnt]['longitude'], $points[$cnt]['latitude']);
830: $cnt++;
831:
832: return $cnt;
833: }
なお、等角航路は複雑な曲線を描くため、距離を求める適当な方程式がない。そこで、隣り合う座標間の大圏距離を積分し、等角航路の距離を同時に算出する。
解説：初期値
greatCircleSailing.php
46: //マップの表示サイズ（単位：ピクセル）
47: define('MAP_WIDTH', 600);
48: define('MAP_HEIGHT', 400);
49: //マップID
50: define('MAPID', 'map_id');
51: //初期値
52: define('DEF_LONGITUDE0', 139.766667); //地図中心（経度）
53: define('DEF_LATITUDE0', 35.681111); // （緯度）
54: define('DEF_FROM', '東京駅'); //from（クエリ）
55: define('DEF_TO', 'ワシントン'); //to （クエリ）
56: define('DEF_CAT_FROM', 'landmark'); //from（カテゴリ）
57: define('DEF_CAT_TO', 'world'); //to （カテゴリ）
58: define('DEF_LONGITUDE1', 139.766667); //from（経度）
59: define('DEF_LATITUDE1', 35.681111); // （緯度）
60: define('DEF_LONGITUDE2', -77.0329165); //to （経度）
61: define('DEF_LATITUDE2', 38.8878665); // （緯度）
62: define('DEF_TYPE', 'roadmap'); //マップタイプ
63: define('DEF_ZOOM', 2); //ズーム
greatCircleSailing.php
193: /**
194: * マーカー設定、フィット機能：Googleマップ用スクリプト
195: * @param float $latitude, $longitude 経路の中央点座標（緯度・経度）
196: * @param bool $drag TRUE：ドラッグで再描画する／FALSE：しない（デフォルト）
197: * @return string JavaScript
198: */
199: function jsDragMarker_Gmap($latitude, $longitude, $drag=FALSE) {
200: $submit = $drag ? 'document.myform.submit();' : '';
201:
202: $js =<<< EOT
203: //イベント設定
204: document.getElementById('fit').addEventListener('click', fitting, false);
205:
206: //マーカー設定
207: marker_A.setDraggable(true);
208: marker_A.addListener('dragend', function() {
209: marker_A.setIcon(icon_A);
210: marker_A.setDraggable(true);
211: latlng = marker_A.getPosition();
212: document.getElementById('longitude1').value = latlng.lng();
213: document.getElementById('latitude1').value = latlng.lat();
214: document.getElementById('from').value = '';
215: {$submit}
216: });
217:
218: marker_B.setDraggable(true);
219: marker_B.addListener('dragend', function() {
220: marker_B.setIcon(icon_B);
221: marker_B.setDraggable(true);
222: latlng = marker_B.getPosition();
223: document.getElementById('longitude2').value = latlng.lng();
224: document.getElementById('latitude2').value = latlng.lat();
225: document.getElementById('to').value = '';
226: {$submit}
227: });
228:
229: /**
230: * 地図を最適サイズにフィットする
231: * @param なし
232: * @return なし
233: */
234: function fitting() {
235: var lat0, lng0, lat1, lng1;
236: var lat = Array();
237: var lng = Array();
238: lat[0] = document.getElementById('latitude1').value;
239: lng[0] = document.getElementById('longitude1').value;
240: lat[1] = document.getElementById('latitude2').value;
241: lng[1] = document.getElementById('longitude2').value;
242: lat[2] = {$latitude};
243: lng[2] = {$longitude};
244:
245: //南西端(lat0, lng0), 北東端(lat1, lng1)を求める
246: if (lat[0] <= lat[1]) {
247: lat0 = lat[0];
248: lat1 = lat[1];
249: } else {
250: lat0 = lat[1];
251: lat1 = lat[0];
252: }
253: if (lat[2] >= 0) {
254: if (lat[0] <= lat[2]) {
255: lng0 = lng[0];
256: lng1 = lng[1];
257: } else {
258: lng0 = lng[1];
259: lng1 = lng[0];
260: }
261: } else {
262: if (lat[0] <= lat[2]) {
263: lng0 = lng[1];
264: lng1 = lng[0];
265: } else {
266: lng0 = lng[0];
267: lng1 = lng[1];
268: }
269: }
270:
271: var bounds = new google.maps.LatLngBounds(new google.maps.LatLng(lat0, lng0), new google.maps.LatLng(lat1, lng1));
272: map.fitBounds(bounds); //地図表示の最適化
273: }
274:
275: EOT;
276: return $js;
277: }
greatCircleSailing.php
429: /**
430: * マーカー設定、フィット機能スクリプト：Leafletマップ用
431: * @param float $latitude, $longitude 経路の中央点座標（緯度・経度）
432: * @param bool $drag TRUE：ドラッグで再描画する／FALSE：しない（デフォルト）
433: * @return string JavaScript
434: */
435: function jsDragMarker_Leaflet($latitude, $longitude, $drag=FALSE) {
436: //Leafletマップ：東経・西経の境目で線分が折り返さないよう調整
437: if ($longitude < 0) $longitude = 360 + $longitude;
438: $submit = $drag ? 'document.myform.submit();' : '';
439:
440: $js =<<< EOT
441: //イベント設定
442: document.getElementById('fit').addEventListener('click', fitting, false);
443:
444: //マーカー設定
445: marker_A.dragging.enable();
446: marker_A.on('dragend', function(event) {
447: var marker = event.target;
448: var position = marker_A.getLatLng();
449: marker.setLatLng(new L.LatLng(position.lat, position.lng),{draggable:true});
450: map.panTo(new L.LatLng(position.lat, position.lng))
451: document.getElementById('longitude1').value = position.lng;
452: document.getElementById('latitude1').value = position.lat;
453: document.getElementById('from').value = '';
454: {$submit}
455: });
456:
457: marker_B.dragging.enable();
458: marker_B.on('dragend', function(event) {
459: var marker = event.target;
460: var position = marker.getLatLng();
461: marker.setLatLng(new L.LatLng(position.lat, position.lng),{draggable:true});
462: map.panTo(new L.LatLng(position.lat, position.lng))
463: document.getElementById('longitude2').value = position.lng;
464: document.getElementById('latitude2').value = position.lat;
465: document.getElementById('to').value = '';
466: {$submit}
467: });
468:
469: /**
470: * 地図を最適サイズにフィットする
471: * @param なし
472: * @return なし
473: */
474: function fitting() {
475: var lat0, lng0, lat1, lng1;
476: var lat = Array();
477: var lng = Array();
478: lat[0] = document.getElementById('latitude1').value;
479: lng[0] = document.getElementById('longitude1').value;
480: lat[1] = document.getElementById('latitude2').value;
481: lng[1] = document.getElementById('longitude2').value;
482: lat[2] = {$latitude};
483: lng[2] = {$longitude};
484:
485: //南西端(lat0, lng0), 北東端(lat1, lng1)を求める
486: if (lat[0] <= lat[1]) {
487: lat0 = lat[0];
488: lat1 = lat[1];
489: } else {
490: lat0 = lat[1];
491: lat1 = lat[0];
492: }
493: if (lat[2] >= 0) {
494: if (lat[0] <= lat[2]) {
495: lng0 = lng[0];
496: lng1 = lng[1];
497: } else {
498: lng0 = lng[1];
499: lng1 = lng[0];
500: }
501: } else {
502: if (lat[0] <= lat[2]) {
503: lng0 = lng[1];
504: lng1 = lng[0];
505: } else {
506: lng0 = lng[0];
507: lng1 = lng[1];
508: }
509: }
510:
511: //Leafletマップ：東経・西経の境目で線分が折り返さないよう調整
512: //bug:Leafletマップではfitが効かない
513: if (lng0 < 0) lng0 = 360 + lng0;
514: if (lng1 < 0) lng1 = 360 + lng1;
515: var bounds = L.latLngBounds([lat0, lng0], [lat1, lng1])
516: map.fitBounds(bounds); //地図表示の最適化
517: }
518: EOT;
519: return $js;
520: }
解説：航路描画スクリプト
greatCircleSailing.php
279: /**
280: * 航路描画スクリプト：Googleマップ用スクリプト
281: * @param int $id 識別ID
282: * @param array $items 航路の座標
283: * @param string $color 描画色（RGB）
284: * @param int $weight 線の太さ（省略時=1）
285: * @return string JavaScript
286: */
287: function jsDrawLine_Gmap($id, $items, $color, $weight=1) {
288: $js =<<< EOT
289: var pline{$id} = new google.maps.Polyline({
290: map: map,
291: path:[
292:
293: EOT;
294: foreach ($items as $item) {
295: $js .=<<< EOT
296: new google.maps.LatLng({$item['latitude']}, {$item['longitude']}),
297:
298: EOT;
299: }
300: $js .=<<< EOT
301: ],
302: strokeColor: "{$color}",
303: strokeOpacity: 1,
304: strokeWeight: {$weight},
305: zIndex: 1
306: });
307:
308: EOT;
309: return $js;
310: }
greatCircleSailing.php
522: /**
523: * 航路描画スクリプト：Leafletマップ用
524: * @param int $id 識別ID
525: * @param array $items 航路の座標
526: * @param string $color 描画色（RGB）
527: * @param int $weight 線の太さ（省略時=1）
528: * @return string JavaScript
529: */
530: function jsDrawLine_Leaflet($id, $items, $color, $weight=1) {
531: $cnt = 0;
532: $js =<<< EOT
533: var pline{$id}_{$cnt} = L.polyline([
534:
535: EOT;
536: foreach ($items as $key=>$item) {
537: if ($key > 0) {
538: //Leafletマップ：東経・西経の境目で線分が折り返さないよう調整
539: if ((($items[$key - 1]['longitude'] > 0) && ($items[$key]['longitude'] < 0)) || (($items[$key - 1]['longitude'] < 0) && ($items[$key]['longitude'] > 0))) {
540: $cnt++;
541: $js .=<<< EOT
542: ], {
543: "color": "#{$color}",
544: "opacity": 1.0,
545: "weight": {$weight},
546: }).addTo(map);
547: var pline{$id}_{$cnt} = L.polyline([
548:
549: EOT;
550: }
551: }
552: $long = ($item['longitude'] < 0) ? 360 + $item['longitude'] : $item['longitude'];
553: $js .=<<< EOT
554: [{$item['latitude']}, {$long}],
555:
556: EOT;
557: }
558: $js .=<<< EOT
559: ], {
560: "color": "#{$color}",
561: "opacity": 1.0,
562: "weight": {$weight},
563: }).addTo(map);
564:
565: EOT;
566: return $js;
567: }
なお、Lefletのメソッド polyline を使って線分を繋げていくと、東経と西経を
跨ぐ地点（経度の符合が変わる地点）を境目に、線分が折り返してしまう。そこで、境目で線分を、いったん断絶させるようにしている。航路を拡大してみると、この部分で切れていることが見えてしまうが、ご了承いただきたい。
また、マーカーの位置を調整するため、前述のユーザー関数 jsDragMarker_Leaflet を含め、マイナス表記の西経を360度表記に変換したりしており、フィット機能が正常に働かないため、これを省いている。
解説：距離計算表
greatCircleSailing.php
570: /**
571: * 距離計算表を作成
572: * @param array $items 地点・距離情報
573: * @param object $pgc pahooGeoCodeクラス
574: * @return string 距離計算表（HTML）
575: */
576: function makeDistanceTable($items, $pgc) {
577: $d11 = sprintf('%.1f', round($items['dist1'], 1));
578: $d12 = sprintf('%.1f', round($pgc->km2mi($items['dist1']), 1));
579: $d13 = sprintf('%.1f', round($pgc->km2nm($items['dist1']), 1));
580: $d21 = sprintf('%.1f', round($items['dist2'], 1));
581: $d22 = sprintf('%.1f', round($pgc->km2mi($items['dist2']), 1));
582: $d23 = sprintf('%.1f', round($pgc->km2nm($items['dist2']), 1));
583: $ww = MAP_WIDTH / 4;
584:
585: $html =<<< EOT
586: <table class="plists">
587: <caption>距離計算</caption>
588: <tr>
589: <th style="width:{$ww}px;"> </th>
590: <th style="width:{$ww}px;">キロメートル</th>
591: <th style="width:{$ww}px;">マイル</th>
592: <th style="width:{$ww}px;">海里</th>
593: </tr>
594: <tr>
595: <th>大圏航路</th>
596: <td>{$d11}</td>
597: <td>{$d12}</td>
598: <td>{$d13}</td>
599: </tr>
600: <tr>
601: <th>等角航路</th>
602: <td>{$d21}</td>
603: <td>{$d22}</td>
604: <td>{$d23}</td>
605: </tr>
606: </table>
607:
608: EOT;
609: return $html;
610: }
活用例
参考サイト
- 各種WebAPIの登録方法：ぱふぅ家のホームページ
- PHPでGoogleを利用して住所から緯度・経度を求める
- PHPで2地点間の直線距離を求める
- 地図上で2点間の直線距離を測：みんなの知識 ちょっと便利帳
（2021年2月20日）PHP8対応，Yahoo! JavaScriptマップ廃止
（2020年3月22日）OSM Nominatim Search APIを利用できるようにした。
（2020年3月14日）地理院地図・OSMに対応した。
（2019年4月25日）フィット機能を追加。IE11動作不具合改善、その他細かい修正を行った。