「地球上の2点間の最短距離」を求めるPHPコード:【ヒュベニの公式】緯度+経度で座標指定。地球儀の球面上で,最短航路の距離[m]+方角[θ]を算出するアルゴリズムのサンプルプログラム

世界地図で2つの座標(緯度と経度)が与えられた時,その2点間の「最短距離」や「方角」を計算して求める方法です。 「XY座標でピタゴラスの定理…?」と考えてはだめですよ! 地球は丸いので,平面幾何は使えないんです。 球面幾何またはヒュベニの公式(Hubenyの式)を使えば,球面上の最短経路の長さや向きがわかります。 PHPのサンプルコードで計算してみましょう。
学問 プログラミング サンプルコード 数学・幾何学・数値計算 アルゴリズム 地理・世界地図 テクノロジー PHP 学術たん 語学たん
11
地球上で「2点間の最短距離とその方角(東西南北)」を,緯度・経度から求める方法は?
ルワンダ語たん(キニャルワンダ語) @rwanda_go_tan
語学たん達の マップ上での位置が、どんどん決まりつつあります。 地図上での位置、つまり 緯度と 経度が決まる。 という事です。 そうなると次は 「この語学たんに、一番近い位置にいる語学たんは誰?」 みたいに ご近所さんを探したくなりますね。 地理的に近い順ランキング を出すわけです。
曲面上で「三平方の定理」は使えない!
ルワンダ語たん(キニャルワンダ語) @rwanda_go_tan
二人の語学たん・方言たんが それぞれ、 異なる緯度・経度をもっています。 この時、二人の間の距離を キロメートル単位で算出する方法は? 平べったい地図なら、ピタゴラスの定理で すぐに距離が出せますが 地球は丸いですからね…。 地球儀という球面、 ないし楕円体の上での距離を計算します。
ルワンダ語たん(キニャルワンダ語) @rwanda_go_tan
2点の座標(緯度・経度) が与えられた時に 2点間の距離の算出は、 大変な数式になります。 Vincentyの式を使えば、 遠く離れた外国間の距離も ミリメートル単位で計算できますが… Vincenty solutions of geodesics on the ellipsoid in JavaScript | Movable Type Scripts movable-type.co.uk/scripts/latlon… pic.twitter.com/5DTDUG53pE
拡大
ルワンダ語たん(キニャルワンダ語) @rwanda_go_tan
もうちょっと簡単な計算式だと、 国土地理院も使っている(であろう) ヒュベニの公式(Hubenyの式) があります。 これは楽に使えますね! 公式の結果だけ使うぶんには、出てくる数学は ラジアンとルートぐらいです。 2つの座標間の距離を求める - Qiita qiita.com/chiyoyo/items/…
「ヒュベニの公式」サンプルコードと動作デモ

地球上(地球儀の球面上)で
実際に最短距離と方角を計算する
PHPのサンプルプログラム
です。

ここでは
「ハワイ語たん」と
「パプア諸語たん」の
球面上での最短距離
を算出しています。

<?php


// ヒュベニの公式で,地球上の2点間の最短距離とその方角を求めるコード
// ※日付変更線をまたぐ場合も正確に動作する
//
// @author:rwanda_go_tan




// ---------- 最短航路を計算する関数 ---------- 


// 2点間の最短距離をヒュベニの公式で求める。(メートル単位)
// 引数の緯度・経度は -180 ~ +180 の範囲内で渡すこと
// @author:rwanda_go_tan
function get_hubeny_distance( 

  $latitude1_deg, $longitude1_deg, // 座標1

  $latitude2_deg, $longitude2_deg  // 座標2

)
{
    // ---------- 緯度・経度の加工 ---------- 

    // 緯度・経度をラジアンに変換
    $latitude1_rad  = deg2rad( $latitude1_deg );
    $longitude1_rad = deg2rad( $longitude1_deg );
    $latitude2_rad  = deg2rad( $latitude2_deg );
    $longitude2_rad = deg2rad( $longitude2_deg );

    // 緯度の差分(ラジアン)
    $diff_latitude_rad = abs( $latitude1_rad - $latitude2_rad );

    // 経度の差算(ラジアン):
    // 日付変更線をまたぐ場合にも
    // 必ず短いほうの弧長を計算するように
    $diff_longitude_rad = min([

      // 単純な差分で済む場合
      abs( $longitude1_rad - $longitude2_rad ),

      // 単純な差分だけでは長くなってしまう場合,
      // 残りの円弧が最短経路となる
      2*pi() - abs( $longitude1_rad - $longitude2_rad )

    ]);

    // 平均緯度(ラジアン)
    $average_latitude_rad = ( $latitude1_rad + $latitude2_rad ) / 2.0;


    // ---------- 各種定数の計算 ---------- 

    // 地球の長半径,赤道半径(Semi-major Axis of the Earth)
    $semi_major_axis = 6378137.0; 

    // 地球の短半径,極半径(Semi-minor Axis of the Earth)
    $semi_minor_axis = 6356752.314245; 

      // 定数の引用元:
      // https://www.trail-note.net/tech/calc_distance/
      // 世界的に使用されている WGS84 (1984年) を採用

    // 離心率を求める
    $E = sqrt(
      ( 
        // 長半径の2乗
        $semi_major_axis * $semi_major_axis 
        - 
        // 短半径の2乗
        $semi_minor_axis * $semi_minor_axis 
      ) 
      / 
      // 長半径の2乗
      ( $semi_major_axis * $semi_major_axis )
    );

    // 子午線(しごせん)曲率半径 M を求める
    $sin_lat = sin( $average_latitude_rad );
    $W = sqrt(
      1.0 
      - 
      $E * $E * ( $sin_lat * $sin_lat )
    );
    $M = $semi_major_axis * ( 1 - $E * $E ) / ( $W * $W * $W );

    // 卯酉線(ぼうゆうせん)曲率半径 N 
    $N = $semi_major_axis / $W; 


    // ---------- 距離算出 ---------- 

    $t1 = $M * $diff_latitude_rad;
    $t2 = $N * cos( $average_latitude_rad ) * $diff_longitude_rad;
    $result = sqrt(
      ( $t1 * $t1 ) 
      + 
      ( $t2 * $t2 )
    );


    // メートル単位で返却
    return $result;
}
// NOTE: Qiitaに載っているサンプルコードは
// 2点の最短航路が日付変更線をまたぐ場合に計算結果がおかしくなる。
// 上記コードではそのバグを修正してある。
// https://qiita.com/chiyoyo/items/b10bd3864f3ce5c56291
// 2点の位置情報から方位角を算出
// 引数の緯度・経度は -180 ~ +180 の範囲で渡す。
// 起点から見た終点へのラジアンを0~2πで返す。北が0で時計回りに角度が増える。
// @author:rwanda_go_tan
function get_azimuth_of_two_points(

  $latitude_from_deg, $longitude_from_deg, // 座標1,起点

  $latitude_to_deg, $longitude_to_deg  // 座標2,終点

){

  // Azimuth : 方位角
  // https://ja.wikipedia.org/wiki/%E6%96%B9%E4%BD%8D%E8%A7%92

  // xの差分
  $diff_longitude_deg = $longitude_to_deg - $longitude_from_deg;

  // NOTE:
  // 緯度・経度は360°単位であり,三角関数での計算時にはラジアンに変換

  // 方位角
  $theta = pi() / 2.0 
    - atan2(

      cos( $latitude_from_deg  / 180.0 * pi() )
      *
      tan( $latitude_to_deg    / 180.0 * pi() )

      -

      sin( $latitude_from_deg  / 180.0 * pi() )
      *
      cos( $diff_longitude_deg / 180.0 * pi() )
      ,

      sin( $diff_longitude_deg / 180.0 * pi() )

    )
  ;

  // 計算式:
  // https://keisan.casio.jp/exec/system/1257670779

  // PHP, atan2( y, x )
  // https://www.php.net/manual/ja/function.atan2.php

  // 角度を0から2πまでに正規化
  while( $theta >= 2*pi() ){ $theta -= 2*pi(); }
  while( $theta < 0       ){ $theta += 2*pi(); }

  return $theta;
}



// 2点の方位角を方角の名称に変換。
// 引数はラジアンで,0~2π。北が0で時計回りに角度が増える。
// @author:rwanda_go_tan
function translate_azimuth_to_direction_name( $azimuth_rad ) {

  // 方角リスト
  $direction_names = [
    "北", "北東", "東", "南東", 
    "南", "南西", "西", "北西", 
    "北" // π/16 だけずらす分を,この末尾要素で回収する
  ];

  // NOTE: 方位を8分割している。
  // つまり,- π/16 ~ + π/16 を「北」とみなす。

  // 方位角を配列のインデックスに変換
  $direction_index = floor(
    (
      // - π/16 を北とみなすために
      ( $azimuth_rad + ( pi() ) / 16.0 )
      /
      ( 2 * pi() )
    )
    *
    8.0 // 2πの回転を,全方位の個数のスキャンに対応させる
  );

  // 方角の名称
  $direction_name = $direction_names[ $direction_index ];

  return $direction_name;
}
// ---------- サンプル座標データ ---------- 


// ハワイ語たん
// https://gogakutan.jishotter.info/detail/hawaii_tan
$hawaigotan = [
  "lat" => 19.8967662,
  "lng" => -155.5827818 
];


// パプア諸語たん
// https://gogakutan.jishotter.info/detail/papuan_tan
$paputan = [
  "lat" => -6.314993,
  "lng" => 143.95555 
];




// ---------- 動作デモ ---------- 



// ヒュベニの公式で最短距離を出す
$distance_hybany = get_hubeny_distance(
  $hawaigotan["lat"], $hawaigotan["lng"], 
  $paputan["lat"],    $paputan["lng"]
);
echo "最短距離 = " . ( $distance_hybany / 1000.0 ) . " km.<br><br><br>";





// 2点の方位角を出す
$azimuth_rad = get_azimuth_of_two_points(
  $hawaigotan["lat"], $hawaigotan["lng"], 
  $paputan["lat"],    $paputan["lng"]
);
echo "方位角 = " . $azimuth_rad . " rad.<br><br><br>";





// 2点の方角名称を出す
$direction_name = translate_azimuth_to_direction_name( $azimuth_rad );
echo "方位名称 = " . $direction_name . " .<br><br><br>";



/*

上記サンプルの出力:

最短距離 = 7285.1919373465 km.
方位角 = 4.4095462906492 rad.
方位名称 = 南西 .

*/


?>
距離計算の用語:大円距離,測地線,方位角
ルワンダ語たん(キニャルワンダ語) @rwanda_go_tan
球面上の最短距離について、 用語の補足… ・地球を球とみなした場合の、球面上の最短距離… 大円距離、大圏距離と呼ばれます。 大円と呼ぶのは、2点を結ぶ最短距離が 2点を含む円周上にあるからですね。 大円距離 ja.m.wikipedia.org/wiki/%E5%A4%A7…
ルワンダ語たん(キニャルワンダ語) @rwanda_go_tan
・(球面上に限らず)曲がった空間内で、2点を結ぶ線… 測地線 ja.m.wikipedia.org/wiki/%E6%B8%AC… 球面上で まっすぐロープを張ったからといって、最短距離とは限りません。 >大円をその上の2点で分けると優弧と劣弧に分かれる。劣弧に沿って移動すれば最短距離、優弧に沿えば直線的な移動としては最も遠回り
ルワンダ語たん(キニャルワンダ語) @rwanda_go_tan
もう1つ知っておきたい用語が 「方位角」… 緯度と経度のわかる 球面上の2点間で、 お互いの存在する方角(東西南北)です。 最短距離を、どの方角に向かって歩けばよいのか 向きを知るためのものですね。 緯度経度から距離と方位を求める logical-arts.jp/archives/18 >方位角θはΔy/Δxの逆正接で求める
ルワンダ語たん(キニャルワンダ語) @rwanda_go_tan
ある語学たんから見て、 別の語学たんはどの方角にいるのか。 角度θがわかれば θを大まかに東西南北に言い換えられますね。 θ=0付近なら東、とか…。 具体的な計算式は 下記などに載っています。 2地点間の距離と方位角 - 高精度計算サイト keisan.casio.jp/exec/system/12… 逆正接=arctanを使います。
地球上の「方角」は,直感とは違った答えになる!
ルワンダ語たん(キニャルワンダ語) @rwanda_go_tan
@Kitakyushutan 角度や方位を 正確に計算するって,面白いテーマで…。 例えば,普通に メルカトル図法の世界地図を見ると 「カナダの イヌクティトゥット語たん @inuktitut_tan から ノルウェーの νのるたん @nynorsk_tan へ 最短距離で向かう方角は "南西" ?」 と思っちゃいますね。 でも答えは 北東なんですpic.twitter.com/wKfFFLlQLM
拡大
ルワンダ語たん(キニャルワンダ語) @rwanda_go_tan
@Kitakyushutan @inuktitut_tan @nynorsk_tan ※南西じゃなくて 南東の打ち間違いでした 地図上ではまっすぐ東か, ノルウェーのほうが やや下側ですもんね。
ルワンダ語たん(キニャルワンダ語) @rwanda_go_tan
@Kitakyushutan @inuktitut_tan @nynorsk_tan 語学たん同士の 方位関係を見てみると… イヌクティトゥット語たん gogakutan.jishotter.info/detail/inuktit…グリーンランド語たん @kalaallisut_tan まで 北東へ1480 ㎞ ・アイスランド語たん @islenska_tan まで 北東へ2787 ㎞ ・νのるたん まで 北東へ4263 ㎞ カナダから,まっすぐ北東へ向かう道なんですね pic.twitter.com/kC0BpRIGMb
拡大
拡大
ルワンダ語たん(キニャルワンダ語) @rwanda_go_tan
@Kitakyushutan @inuktitut_tan @nynorsk_tan @kalaallisut_tan @islenska_tan 本家Google Mapは,地球儀のように 丸い地球を見ることが可能になってます。 北極や南極周辺の リアルな地形は メルカトル図法では把握できず 地球儀で見ると,よくわかりますね! カナダ→グリーンランド→アイスランド→ノルウェー これは,カナダから見ると 北東向きの 一直線の航路なのです。 pic.twitter.com/ZThRbCN8wp
拡大
「ヒュベニの式」世間のサンプルコードには,バグが潜んでいる場合もある!
ルワンダ語たん(キニャルワンダ語) @rwanda_go_tan
@Kitakyushutan @inuktitut_tan @nynorsk_tan @kalaallisut_tan @islenska_tan 地球の球面上で, 2座標の最短距離を求める「ヒュベニの公式」。 実は,世間に出回っている 「ヒュベニの公式」のサンプルコードには, バグがあります。 日付変更線をまたぐ】時に, 正しく距離を算出できないのです! 経度 +180° 側から 経度 -180° 側へまたぐ時とかね。
ルワンダ語たん(キニャルワンダ語) @rwanda_go_tan
@Kitakyushutan @inuktitut_tan @nynorsk_tan @kalaallisut_tan @islenska_tan 「ヒュベニの公式」実装バグの例として ハワイにいる「ハワイ語たん@hawaii_tan と パプアニューギニアにいる 「パプア諸語たん@papuan_tan の 2点間の最短距離を求めてみます なんと,下図・青線のように 「日付変更線をまたがない」長~い距離を 誤って計算してしまう場合があるのです! pic.twitter.com/sFvl3lLbMC
拡大
ルワンダ語たん(キニャルワンダ語) @rwanda_go_tan
@Kitakyushutan @inuktitut_tan @nynorsk_tan @kalaallisut_tan @islenska_tan @hawaii_tan @papuan_tan ハワイ⇔パプアニューギニア間ならば 本当は,下の図の オレンジ色の線が「最短距離」ですよね。 この航路は,日付変更線をまたぎます。 最短距離を計算する際, 流通しているサンプルコードを見ると どうも,その部分の処理の実装が 漏れているケースがあるようなのです。 pic.twitter.com/7lmFMIomRo
拡大
正しく計算できるように,バグを除去しました。
ルワンダ語たん(キニャルワンダ語) @rwanda_go_tan
@Kitakyushutan @inuktitut_tan @nynorsk_tan @kalaallisut_tan @islenska_tan @hawaii_tan @papuan_tan そこで,私は 「日付変更線をまたぐ場合も 正しく最短距離を出せるプログラム」 を作りました。 その結果, 世界中の「語学たん」どうしの距離を 正しく計算できています。 サンプルとして, ハワイ語たんと パプア諸語たんの間の最短距離は 下記ページから見れます: gogakutan.jishotter.info/distance?from=… pic.twitter.com/4Elt2NhCrW
拡大
拡大
残りを読む(20)

コメント

ルワンダ語たん(キニャルワンダ語) @rwanda_go_tan 2019年8月17日
デコレーション機能を駆使して,コード部分をシンタックスハイライトしようと試みましたが,途中で挫折しました。
さろげーと @surrogatepair 2019年8月17日
QiitaのコードにバグがあるのであればQiita側で編集リクエスト送ったほうがよさそう
ログインして広告を非表示にする
ログインして広告を非表示にする