どうもこんにちは。iOSをメインに開発しているロッキーカナイです。
私は現在、地図系アプリを開発しているのですが、緯度経度から距離を求める必要があって色々調べたのですが、Swiftで紹介しているサイトが無かったので、ここで紹介したいと思います。
計算方法
今回は3つの計算方法を紹介します。
- 球面三角法
- ヒュベニの公式
- 測地線航海算法
球面三角法
GoogleMapAPIのgeometory.computeDistanceBetweenでも使用されています。誤差はあるものの計算式が簡素なのが特徴です。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
func distance(current: (la: Double, lo: Double), target: (la: Double, lo: Double)) -> Double { // 緯度経度をラジアンに変換 let currentLa = current.la * Double.pi / 180 let currentLo = current.lo * Double.pi / 180 let targetLa = target.la * Double.pi / 180 let targetLo = target.lo * Double.pi / 180 // 赤道半径 let equatorRadius = 6378137.0; // 算出 let averageLat = (currentLa - targetLa) / 2 let averageLon = (currentLo - targetLo) / 2 let distance = equatorRadius * 2 * asin(sqrt(pow(sin(averageLat), 2) + cos(currentLa) * cos(targetLa) * pow(sin(averageLon), 2))) return distance / 1000 } |
ヒュベニの公式
一般的に使用されている方法で、国土交通省の距離計算にも使われているようです。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 |
func distance2(current: (la: Double, lo: Double), target: (la: Double, lo: Double)) -> Double { // 緯度経度をラジアンに変換 let currentLa = current.la * Double.pi / 180 let currentLo = current.lo * Double.pi / 180 let targetLa = target.la * Double.pi / 180 let targetLo = target.lo * Double.pi / 180 // 緯度差 let radLatDiff = currentLa - targetLa // 経度差算 let radLonDiff = currentLo - targetLo // 平均緯度 let radLatAve = (currentLa + targetLa) / 2.0 // 測地系による値の違い // 赤道半径 // let a = 6378137.0 world let a = 6377397.155 // japan // 極半径 // let b = 6356752.314140356 world let b = 6356078.963 // japan // 第一離心率^2 let e2 = (a * a - b * b) / (a * a) // 赤道上の子午線曲率半径 let a1e2 = a * (1 - e2) let sinLat = sin(radLatAve); let w2 = 1.0 - e2 * (sinLat * sinLat); // 子午線曲率半径m let m = a1e2 / (sqrt(w2) * w2); // 卯酉線曲率半径 n let n = a / sqrt(w2) // 算出 let t1 = m * radLatDiff let t2 = n * cos(radLatAve) * radLonDiff let distance = sqrt((t1 * t1) + (t2 * t2)) return distance / 1000 } |
測地線航海算法
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 |
func distance3(current: (la: Double, lo: Double), target: (la: Double, lo: Double)) -> Double { // 緯度経度をラジアンに変換 let currentLa = current.la * Double.pi / 180 let currentLo = current.lo * Double.pi / 180 let targetLa = target.la * Double.pi / 180 let targetLo = target.lo * Double.pi / 180 // 赤道半径 let a = 6378137.0 // 極半径 let b = 6356752.314140356 // 扁平率 let f = (a - b) / a // 算出 let p1 = atan(b / a) * tan(currentLa) let p2 = atan(b / a) * tan(targetLa) let sd = acos(sin(p1) * sin(p2) + cos(p1) * cos(p2) * cos(currentLo - targetLo)) let cos_sd = cos(sd / 2) let sin_sd = sin(sd / 2) let c = (sin(sd) - sd) * pow(sin(p1) + sin(p2), 2) / cos_sd / cos_sd let s = (sin(sd) + sd) * pow(sin(p1) - sin(p2), 2) / sin_sd / sin_sd let delta = f / 8.0 * (c - s) let distance = a * (sd + delta) return distance / 1000 } |
まとめ
上記の各計算式を用いて比較してみます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
// 現在地の緯度経度 let currentPosition = (la: 35.68944, lo: 139.69167) // ターゲットの緯度経度 let targetPosition = (la: 35.85694, lo: 139.64889) let di = distance(current: currentPosition, target: targetPosition) print("球面三角法 距離:\(di)km") let di2 = distance2(current: currentPosition, target: targetPosition) print("ヒュベニの公式 距離:\(di2)km") let di3 = distance3(current: currentPosition, target: targetPosition) print("測地線航海算法 距離:\(di3)km") |
現在地の緯度経度を東京都庁(緯度:35.68944 経度:139.69167)にし、目的地の緯度経度を埼玉県庁(緯度:35.85694 経度:139.64889)にしました。
県庁所在地の位置情報はこちらより拝借致しました。
各計算式の結果が以下の通りです。
1 2 3 |
球面三角法 距離:19.0421298714453km ヒュベニの公式 距離:18.9811962944882km 測地線航海算法 距離:22.5088536129471km |
ちなみに国土交通省の県庁所在地間の距離のデータがありますので答え合わせをしてみたいと思います。
東京都庁と埼玉県庁の距離は19.0kmとあるので、球面三角法とヒュベニの公式が近い値となりました。もしかしたら測地線航海算法の計算に誤りがあるかも。。。
今開発しているアプリでは計算方法がシンプルな球面三角法を使用する事になりそうです。
ではー
ABOUT ME
