もくじ
どうもこんにちは。iOSをメインに開発しているロッキーカナイです。
私は現在、地図系アプリを開発しているのですが、緯度経度から距離を求める必要があって色々調べたのですが、Swiftで紹介しているサイトが無かったので、ここで紹介したいと思います。
計算方法
今回は３つの計算方法を紹介します。
- 球面三角法
- ヒュベニの公式
- 測地線航海算法
球面三角法
GoogleMapAPIのgeometory.computeDistanceBetweenでも使用されています。誤差はあるものの計算式が簡素なのが特徴です。
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
}ヒュベニの公式
一般的に使用されている方法で、国土交通省の距離計算にも使われているようです。
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
}測地線航海算法
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
}まとめ
上記の各計算式を用いて比較してみます。
// 現在地の緯度経度
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）にしました。
県庁所在地の位置情報はこちらより拝借致しました。
各計算式の結果が以下の通りです。
球面三角法 距離:19.0421298714453km
ヒュベニの公式 距離:18.9811962944882km
測地線航海算法 距離:22.5088536129471kmちなみに国土交通省の県庁所在地間の距離のデータがありますので答え合わせをしてみたいと思います。
東京都庁と埼玉県庁の距離は19.0kmとあるので、球面三角法とヒュベニの公式が近い値となりました。もしかしたら測地線航海算法の計算に誤りがあるかも。。。
今開発しているアプリでは計算方法がシンプルな球面三角法を使用する事になりそうです。
ではー
