[지리 정보] Vincenty 공식을 이용한 지도 좌표 거리 계산 및 Kotlin 함수 작성

- 이번 포스팅으론, 이전 글에서 구현 및 테스트를 진행한 적이 있는 Harversine 공식을 사용한 좌표 거리 계산 함수를 개선하기 위하여 Vincenty 공식을 알아보고 구현하겠습니다.

 

- 앞선 글에서 Harversine 공식은 지구를 완전한 구형 물체로 가정하여 계산을 진행한다는 단점이 있다고 설명드렸습니다.

실제 지구는 완전한 구형이 아닌, 자전 방향으로 조금 튀어나와 있는 미묘한 타원체라고 할 수 있으므로, 해당 방식은 당연히 오차를 만들 수 밖에 없는데,

이번에는 그러한 단점을 보안하여 조금 더 정확도를 높인 Vincenty 공식에 대해 공부해봅시다.

 

- Vincenty 공식의 주요 개념

vincenty 공식은 지구를 타원형으로 본다고 설명드렸죠? 그렇다면 타원을 이루는 두가지 변수가 존재함을 알아야 합니다.

구를 타원체로 만들어주는 두 변수는, 장반경과 단반경입니다.

장반경은 타원체의 긴 축의 반지름이고, 단반경은 짧은 축의 반지름입니다.

이 두 변수에 따라 타원체의 전체 형태가 정해지며,

아래에 설명드릴 Vincenty 공식에서 지구 형태를 가정할 때에는, 좌표계 시스템 중 WGS-84 타원체를 가정하며,

이는 장반경(a)이 6378137.0m(지구 적도 반경), 단반경(b)이 6356752.314245m(지구의 극 반경) 인 타원체입니다.

 

위 a, b 변수를 통해 편평률이라는 값을 계산할 수 있는데,

편평률은 타원체의 장반경과 단반경 간의 차이를 나타내는 값으로, 타원체가 얼마나 납작해져 있는지(편평해진 정도)를 표현합니다. 즉, 완전한 구에서 벗어나 타원이 된 정도를 나타내는 척도입니다.

 

편평률의 계산법은,

f = (a-b) / a

이며, WGS-84 타원체의 편평률은 계산하자면,

(6378137 - 6356752.314245 ) / 6378137 = 0.0033528106647756

로 계산이 되네요.

 

위와 같은 상수값을 사용하여 좌표 계산의 정확도를 높인 vincenty 알고리즘을 구현한 코드는,

    // (지도 좌표 1 에서 지도 좌표 2 까지의 거리 (미터) 반환, Vincenty 공식)
    // Vincenty 공식은 타원체 위에서 두 좌표 사이의 거리를 계산하는 방법으로, 지구를 완전한 구가 아닌 타원체로 간주하기 때문에 더 정확한 결과를 제공합니다.
    // 좌표계는 WGS84 를 사용합니다.
    fun getDistanceMeterBetweenTwoLatLngCoordinateVincenty(
        latlng1: Pair<Double, Double>,
        latlng2: Pair<Double, Double>
    ): Double {
        val a = 6378137.0 // WGS-84 기준 타원체의 장반경 (미터)
        val b = 6356752.314245 // 단반경
        val f = 0.0033528106647756 // WGS-84 타원체의 편평률

        val lat1 = latlng1.first * PI / 180
        val lat2 = latlng2.first * PI / 180
        val lon1 = latlng1.second * PI / 180
        val lon2 = latlng2.second * PI / 180

        val U1 = atan((1 - f) * tan(lat1))
        val U2 = atan((1 - f) * tan(lat2))
        val L = lon2 - lon1
        var lambda = L

        val maxIterations = 200
        val tolerance = 1e-12
        var cosSqAlpha: Double
        var sinSigma: Double
        var cos2SigmaM: Double
        var cosSigma: Double
        var sigma: Double
        var lambdaP: Double
        var iter = 0

        do {
            val sinLambda = sin(lambda)
            val cosLambda = cos(lambda)
            sinSigma = sqrt(
                (cos(U2) * sinLambda) * (cos(U2) * sinLambda) +
                        (cos(U1) * sin(U2) - sin(U1) * cos(U2) * cosLambda) *
                        (cos(U1) * sin(U2) - sin(U1) * cos(U2) * cosLambda)
            )

            if (sinSigma == 0.0) return 0.0 // 두 좌표가 동일함

            cosSigma = sin(U1) * sin(U2) + cos(U1) * cos(U2) * cosLambda
            sigma = atan2(sinSigma, cosSigma)
            val sinAlpha = cos(U1) * cos(U2) * sinLambda / sinSigma
            cosSqAlpha = 1 - sinAlpha * sinAlpha
            cos2SigmaM = if (cosSqAlpha != 0.0) cosSigma - 2 * sin(U1) * sin(U2) / cosSqAlpha else 0.0
            val C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
            lambdaP = lambda
            lambda = L + (1 - C) * f * sinAlpha * (sigma + C * sinSigma * (
                    cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)))
        } while (abs(lambda - lambdaP) > tolerance && ++iter < maxIterations)

        if (iter >= maxIterations) return Double.NaN // 수렴 실패

        val uSq = cosSqAlpha * (a * a - b * b) / (b * b)
        val A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
        val B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
        val deltaSigma = B * sinSigma * (
                cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) -
                        B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) *
                        (-3 + 4 * cos2SigmaM * cos2SigmaM)))

        val s = b * A * (sigma - deltaSigma)

        return s // 거리 (미터)
    }

 

위와 같습니다.

 

코드가 아닌 수식에 대한 해석은 구면 삼각함수 설명, Harversine 수식 설명 이후 정리하도록 하겠습니다.

 

- 이상입니다.