하버사인 공식

출처: 나노 바나나

정의에 따르면 하버사인 공식(Haversine, Half Versed sine, 삼각함수의 일종)은 지구와 같은 구체 위에서 두 지점의 위도와 경도를 이용해 가장 짧은 거리(대원 거리)를 정확하게 계산하는 공식이다. 구면 코사인 법칙을 활용할 수도 있는데 이 공식에서는 두 점의 거리가 매우 가까울 때 정밀도가 떨어지는 단점 때문에 하버사인 공식으로 변형해서 사용한다고 한다. 앞으로 이어질 설명은 수학적 증명이나 유도와 거리가 멀다. 공식을 추상적으로 이해한 뒤 코드 변환에 중점을 둔다.

 

개념적인 부분을 중점적으로 짚어보자. 그 전에 확실하게 해야 할 기본 개념으로 라디안(호도법)이 있다. 일상 생활에서 쓰이는 각도법은 직관적으로 계산하기 편하지만, 그 자체로 호의 길이와 대응되지 않는다. 호도법에서, 반지름이 r일 때 호의 길이가 r을 이루는 각도를 1라디안으로 본다.

 

호의 길이가 2r라면 중심 각도는 2라디안, 호의 길이가 3r이라면 3라디안이 된다. 즉 중심각도가 theta일 때 호의 길이는 r * theta인 셈. 원의 둘레(원주)가 2  * PI * r일 때 중심각도는 360도가 된다. 여기에서 2 * PI가 360도임을 알 수 있다. 1도는 PI/180이 되기 때문에 Degree를 라디안으로 변환할 때, 원하는 각도에 PI/180을 곱하면 된다.

 

처음에 언급했듯 하버사인 공식은 두 위경도 좌표 사이의 거리를 계산할 때 사용한다. 지도 위의 도로를 이동할 때 최단 거리를 구하려면 간선의 가중치를 알아야 한다. 하버사인 공식을 통해 두 지점 간에 거리를 구해 가중치로 사용한다.

 

지구와 같은 구면 위에 그려진 삼각형(구면 삼각형)의 변과 각도 사이의 관계를 나타내는 기본 법칙으로 구면 코사인 법칙이 있다. 구면 위의 세 점 A, B, C가 있고 그 대변의 길이(각도)를 a, b, c라고 할 때 다음 공식이 성립한다.

 

여기서 a, b, c가 길이(각도)라는 표현이 나오는데 "길이 = 각도"는 직관적으로 이해하기 어렵다(이미 알고 있는 내용이면, 이 글을 읽을 필요 없을지도). a만 놓고 보자. 점 A의 대변(실제로는 호)이고 이에 대한 중심각은 지구 중심과 이루는 각도다. 그리고 이 모든 계산 과정에서 반지름이 1인 단위 원을 전제로 해야 한다.

 

라디안 설명에서, 1라디안 일 때 호의 길이를 r, 2라디안일 때 2r, 3라디안일 때 3r이라고 했었다. 호의 길이는 theta x r(반지름)이라는 공식을 얻을 수 있고 단위 원에서는 반지름이 1이기 때문에 길이 = 각도의 공식이 성립된다.

 

구면 코사인 법칙을 사용하면 두 지점이 너무 가까울 때 컴퓨터에서 부동소수점을 계산하는 방식 때문에 정밀도가 떨어지는 단점이 있다. 이를 보완하기 위해 하버사인 공식을 사용할 수 있다. 하버사인 공식으로 넘어가기 전에 위의 구면 코사인 법칙을 지구 상의 두 지점(우리가 궁극적으로 구해야 하는 거리)과 북극점을 잇는 구면 삼각형을 보자.

b의 길이는 위도 90⁰ - 위도A가 되며 a의 길이는 위도 90⁰ - 위도B가 된다. 이 값을 구면 코사인 법칙에 넣을 수 있는데 수학에서 cos(90⁰ - θ) = sin(θ)이고 sin(90 - θ) = cos(θ) 라는 성질이 있다. 그리고 끼인각 C(∠ANB)는 A와 B의 경도 차이인 ∆λ(람다)로 표현하면 다음으로 정리된다.

다음부터는 그냥 간단하게 넘어갈 건데(must), 반각 공식인 cos(θ) = 1 - 2sin²(θ/2)을 적용해서 다음 공식을 유도할 수 있다고 한다.

 

길을 잃으면 안 된다. 우리는 c를 구하는 것이 목적이다. 우변은 우리가 적정한 값을 넣어서 값으로 구할 수 있다. 그렇게 구한 값에 루트(√)를 씌우고 아크사인을 취하면 이론적으로 중심각 c를 구할 수 있다.

아크사인을 쓰지 않을 건데 그 이유는 두 지점 A, B가 만약 지구 정반대편에 있어 거리가 가장 멀 때, 부동소수점 연산이 한계로 1이 아니라 1.00000000000002 같은 값이 나오면 Math.asin() 함수가 NaN을 반환한다. 아크탄젠트를 이용한다.

 

아크탄젠트는 밑변과 높이를 통해 각도를 알려준다. 직각 삼각형을 만들기 위해 빗변의 길이가 1(반지름 1)이고 각도가 c/2인 직각 삼각형을 만들어보자. sin(c/2)가 높이가 되고 위의 수식으로 sin(c/2)는 √a인 걸 알고 있다. 밑변은 피타고라스의 정의로 구할 수 있다.

 

그렇게 해서 아크탄젠트에 높이와 밑변을 넣어서 구한 값은 c/2 기준이기 때문에 2를 곱하면 우리가 원하는 각도(길이)가 나온다. 참고로 관련 함수는 Math.atan()이 아니라 Math.atan2()를 사용한다. 이유는, Math.atan()은 나눗셈 연산이 들어가므로 0으로 나눌 위험이 있고 Math.atan2()는 내부적으로 0이 되는 상황을 안전하게 예외처리하도록 설계되어 있다. 우리가 구한 각도(길이)는 단위 원 기준이기 때문에 지구의 반지름을 곱해 정확한 거리를 구할 수 있다. 다만 우리는 이 거리 값을 가중치로 사용할 것이므로 적당한 값을 곱해도 된다.

 

이 모든 과정을 코드로 표현하면 다음과 같다.

function getDistanceFromLatLonInM(lon1, lat1, lon2, lat2) {
  const R = 6371; // 지구의 반지름 (킬로미터)
  const dLat = (lat2 - lat1) * (Math.PI / 180);
  const dLon = (lon2 - lon1) * (Math.PI / 180);

  // sin^2(c/2)
  const a =
    Math.sin(dLat / 2) * Math.sin(dLat / 2) +
    // 경선의 간격이 극으로 가면서 좁아지는 걸 반영하도록 보정
    Math.cos(lat1 * (Math.PI / 180)) *
      Math.cos(lat2 * (Math.PI / 180)) *
      //
      Math.sin(dLon / 2) *
      Math.sin(dLon / 2);

  // sin^2(theta) + cos^2(theta) = 1의 법칙을 활용  
  const c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));

  // 반지름에 중심각을 곱하면 호의 길이
  return R * c;
}