【Kotlin】Hubenyの公式で距離計算
Last-modified: 2020-04-09 (木) 20:49:06 (1439d)
Hubenyの公式で距離計算 †
2点の地理座標間の距離を計算するのに、そんなに複雑ではなくて、球面三角法よりももう少し精度の良い計算でHubenyの公式というものがあるみたいです。
とりあえず簡単にKotolinで試してみます。
精度は、日本国内程度の近距離ならそんなに悪くなさそう(1%以下)だけれど、赤道をまたいだり半球の1/4を超えるような距離だと結構ずれる様子(3-5%くらい?)。大圏コースとか測地線とか関係あるのかもしれない。
// https://ja.wikipedia.org/wiki/%E6%B8%AC%E5%9C%B0%E5%AD%A6 val WGS84_A = 6378137.0 //長半径 val WGS84_B = 6356752.314245 //短半径 val WGS84_E2 = ((WGS84_A * WGS84_A) - (WGS84_B * WGS84_B)) / (WGS84_A * WGS84_A) // 離心率 val WGS84_A1E2 = WGS84_A * (1 - WGS84_E2) // 赤道上の子午線曲率半径 fun deg2rad(d: Double): Double { return d * Math.PI / 180.0 } // 小数点で緯度経度を取り、距離(m)を返す。 fun calcHubeny(lat1: Double, lng1: Double, lat2: Double, lng2: Double): Double { val latdiff = deg2rad(lat1 - lat2) // 緯度差 val lngdiff = deg2rad(lng1 - lng2) // 経度差 val latavg = deg2rad((lat1 + lat2) / 2.0) // 緯度平均 //卯酉線曲率半径 val sinLatAvg = Math.sin(latavg) val w2 = 1.0 - WGS84_E2 * (sinLatAvg * sinLatAvg) val n = WGS84_A / Math.sqrt(w2) //子午線曲率半径 val m = WGS84_A1E2 / (Math.sqrt(w2) * w2) // Hubeny val t1 = m * latdiff val t2 = n * Math.cos(latavg) * lngdiff return Math.sqrt((t1 * t1) + (t2 * t2)) }
2点の座標から方角を求める †
/*** * 小数点で2点の緯度経度から、1から見た2の方角を, * 北を0度とした360度の整数で返す */ fun getDirection (lat1: Double, lng1: Double, lat2: Double, lng2: Double):Int{ val Y = Math.cos(lng2 * Math.PI / 180.0) * Math.sin(lat2 * Math.PI / 180.0 - lat1 * Math.PI / 180.0) val X = Math.cos(lng1 * Math.PI / 180.0) * Math.sin(lng2 * Math.PI / 180.0) - Math.sin(lng1 * Math.PI / 180.0) * Math.cos(lng2 * Math.PI / 180.0) * Math.cos(lat2 * Math.PI / 180.0 - lat1 * Math.PI / 180.0) var E0 = 180.0 * Math.atan2(Y, X) / Math.PI if(E0 < 0){ E0 += 360.0 } val N0 = ((E0 + 90.0) % 360).roundToInt() return N0 }
おまけ †
android用
<string name="n">北</string> <string name="nne">北北東</string> <string name="ne">北東</string> <string name="ene">東北東</string> <string name="e">東</string> <string name="ese">東南東</string> <string name="se">南東</string> <string name="sse">南南東</string> <string name="s">南</string> <string name="ssw">南南西</string> <string name="sw">南西</string> <string name="wsw">西南西</string> <string name="w">西</string> <string name="wnw">西北西</string> <string name="nw">北西</string> <string name="nnw">北北西</string>
val DIREC = arrayOf( getString(R.string.n), getString(R.string.nne), getString(R.string.ne), getString(R.string.ene), getString(R.string.e), getString(R.string.ese), getString(R.string.se), getString(R.string.sse), getString(R.string.s), getString(R.string.ssw), getString(R.string.sw), getString(R.string.wsw), getString(R.string.w), getString(R.string.wnw), getString(R.string.nw), getString(R.string.nnw) ) fun get16Direction(d: Double):String { val i = (((d + 11.25) % 360)/22.5).toInt() return DIREC[i] }