[VB.NET] [Algorithm] 지구에서 두 곳의 위치 사이의 거리 구하기

지구에서 두 점 사이의 거리를 구하는 방법을 알아보도록 합시다.

Haversine (하버사인) 공식을 이용하여 최단거리 구하기.

Haversine (하버사인) 공식은 구에서 두 점 사이의 거리를 계산할 때 사용되는 공식입니다.
다시 말해, 지구 표면에서 가장 짧은 거리인 위치들 사이의 합을 도출해낼 수 있습니다.
(당연하겠지만, 지구의 형태가 완벽한 구 형태가 아니므로 미세한 오차가 존재하긴 합니다.)

아래 예시들은 모두 비주얼베이직 (VB.NET) 언어로 작성되어 있습니다만, 대부분의 프로그래밍 언어에서 제공하는 Math 함수를 이용한 식과 크게 다르지 않습니다.

Haversine Formula

두 위도, 경도 좌표를 기준으로 두 지점 사이의 거리를 구할 때 사용하는 것이 하버사인 공식입니다. 공식은 다음과 같습니다.

θ 는 두 점을 잇는 호의 중심각입니다. (단위는 Radian 입니다)

φ1, φ2 : 1 지점과 2 지점의 위도 (단위는 Radian 입니다)
λ1, λ2 : 1 지점과 2 지점의 경도 (단위는 Radian 입니다)

hav(θ) 는 Haversine (하버사인) 함수로 다음과 같이 표현됩니다.

거리를 구하기 위해서 역함수인 Arc-Haversine (아크-하버사인) 을 곱해줍니다.

즉, 마지막 공식을 이용하면 두 지점 사이의 거리를 구할 수 있습니다.

a = sin²(Δφ / 2) + cos φ1 × cos φ2 × sin²(Δλ / 2)
c = 2 × atan2(√a, √(1−a))
d = R × c

φ 는 latitude (위도), λ 는 longitude (경도), R 은 지구의 반지름 (6,371 km) 입니다.

소스 코드

            Dim pi1 = latx * Math.PI / 180
            Dim pi2 = laty * Math.PI / 180
            Dim triPi = (laty - latx) * Math.PI / 180
            Dim DeltaLambda = (lngy - lngx) * Math.PI / 180

            Dim x As Double
            Dim y As Double
            Dim z As Double

            x = Math.Sin(triPi / 2) * Math.Sin(triPi / 2) + Math.Cos(pi1) * Math.Cos(pi2) * Math.Sin(DeltaLambda / 2) * Math.Sin(DeltaLambda / 2)
            y = 2 * Math.Atan2(Math.Sqrt(x), Math.Sqrt(1 - x))
            z = r * y

            'MessageBox.Show(Val(z / 1000))
            MessageBox.Show(Math.Round(Val(z / 1000), 5) & " km")

Haversine (하버사인) 공식은 비교적 짧은 거리에서도 수치 계산을 위한 양호한 상태를 유지하므로, Spherical Law of Cosines (구면 코사인 법칙) 에 기초한 계산과는 계산 과정과 결과 값이 다릅니다. a 는 두 점 사이의 현 길이 절반의 제곱이며, c 는 각각의 거리를 Radian 으로 표시한 값입니다.

만약 atan2 함수를 사용할 수 없으면, c 를 2 × asin( min(1, √a) ) 로 계산할 수도 있습니다.
(단, 반올림 시 오류 발생에 대한 별도의 처리가 필요합니다.)

Spherical Law of Cosines (구면 코사인 법칙) 을 이용하여 최단거리 구하기.

대부분의 최신 PC 및 프로그래밍 언어에서는 ‘IEEE 754’ 규격을 지원하는 관계로 64-bit의 부동 소수점 수를 사용하며, 소수점 아래 최대 15 자리까지 표현할 수 있습니다.

이에 따른 가설에 의하면, 신뢰성이 높은 정확도로 코사인 공식 (cos c = cos a × cos b + sin a × sin b × cos c) 의 간단한 구면 법칙은 지구 표면에서 몇 m 만큼의 작은 거리까지 잘 보정된 결과를 제공할 수 있습니다. 다만, 상황에 따라서는 아주 작은 거리의 경우에 한해 정방형 근사법을 적용하는 방법이 더 적합 할 수 있습니다.

이렇게, 코사인의 간단한 법칙을 통해 대량의 측지 목적에 한해, Haversine (하버사인) 공식에 대한 합리적인 대안으로 볼 수 있습니다.

Law of Cosines

d = acos(sin φ1 × sin φ2 + cos φ1 × cos φ2 × cos Δλ) × R

소스 코드

            Dim pi1 = latx * Math.PI / 180
            Dim pi2 = laty * Math.PI / 180
            Dim DeltaLambda = (lngy - lngx) * Math.PI / 180

            Dim z = Math.Acos(Math.Sin(pi1) * Math.Sin(pi2) +
                    Math.Cos(pi1) * Math.Cos(pi2) * Math.Cos(DeltaLambda)) * r

            'MessageBox.Show(Val(z / 1000))
            MessageBox.Show(Math.Round(Val(z / 1000), 5) & " km")

Haversine (하버사인) 공식보다 계산과정이 간단한 반면, 비교적으로 계산속도는 약간 느릴 수 있습니다.

Equirectangular Approximation (정방형 근사)

계산 성능이 조금 더 높고, 신뢰성이 중요하지 않은 상황이라면, 피타고라스 정리를 정방형 투영에 활용하실 수 있습니다.

Equirectangular Approximation

x = Δλ × cos φ𝚖
y = Δφ
d = R × √x² + y²

소스 코드

            Dim pi1 = latx * Math.PI / 180
            Dim pi2 = laty * Math.PI / 180
            Dim DeltaLambda = (lngy - lngx) * Math.PI / 180

            Dim x = DeltaLambda * Math.Cos((pi1 + pi2) / 2)
            Dim y = (pi2 - pi1)
            Dim z = Math.Sqrt(x * x + y * y) * r

            'MessageBox.Show(Val(z / 1000))
            MessageBox.Show(Math.Round(Val(z / 1000), 5) & " km")

이 방법으로는 한 번의 삼각함수 (cosine) 와 1 회의 제곱근 (square root) 으로 값을 도출해내는 방법입니다.
앞서, Haversine (하버사인) 공식이 7 회의 삼각함수와 2 회의 제곱근 (square root) 을 사용하는 것과는 차이가 있습니다.

정확도나 신뢰성은 상대적으로 떨어지며, 자오선 (지평면의 남북점, 천정, 천저를 지나는 선) 을 따른다면 계산 오차가 크지 않습니다. 그렇지 않으면 거리, 방위 위도에 따라 계산결과가 달라지지만 사소한 수준입니다. 앞서 말씀드렸듯이, 좌표간의 거리가 길수록 오차율은 상당히 높아집니다.