위도, 경도 값을 가지고 두 지점의 거리를 구하는 일은 일반적으로 접하지 않을 듯 하지만, 기록 차원에서 남겨 놓습니다. 아마 지도 api를 이용해서 프로그램 등을 개발을 하시는 분들이 아니라면 언제 써먹을 일이 있을까 싶습니다.
개인적으로는 업무상 특정 이벤트가 발생한 지점과 특정 시설 사이의 거리를 계산해야 할 일이 있었습니다. 구글맵이나 네이버/카카오 지도를 웹으로 열어 거리를 측정할 수 있긴 합니다. 그것을 하나하나 수작업으로 할 것이 아니라, 위도 및 경도 값을 받아서 여러 장소에 대해 계산을 해야 했습니다.
위도 및 경도 값을 이용해서 두 지점 사이의 거리를 구하는 방법은 흔히 하버사인 공식(Haversine Formula)이 주로 나옵니다. 지구를 완벽한 구형으로 가정하여 계산을 하는 것이라 매우 단순합니다. 처음에는 저도 하버사인 공식을 이용했었으나, 거리가 멀어질수록 값이 꽤 차이가 있었습니다. 엄청 정밀한 것을 원한 것은 아니었으나, 공식적인 값하고 경우에 따라 수 km씩 오차가 생기기도 했습니다.
이래저래 구글링을 해보니 WGS-84 회전 타원체라는 개념도 알게 되고, Vincenty formula를 알게 되었습니다. 이들 개념과 이를 이용한 Golang 코드를 공유합니다.
WGS-84 타원체
- 정의: WGS-84(World Geodetic System 1984) 타원체는 지구를 근사하는 회전 타원체 모델입니다. 이 모델은 지구의 평균적인 모양과 크기를 나타내며, GPS와 기타 위성 기반 위치 서비스에서 표준으로 사용됩니다.
- 특징: WGS-84 타원체는 두 주요 매개변수로 정의됩니다:
- 장반경(a): 타원체의 중심에서 적도까지의 거리(약 6,378.137 킬로미터).
- 편평률(f): 타원체의 압축을 나타내며, 장반경과 단반경(b) 사이의 비율로 표현됩니다. WGS-84의 편평률은 약 1/298.257223563입니다.
- 용도: 이 타원체는 지구의 표면을 모델링하고, 위도와 경도를 사용하여 지구상의 위치를 결정하는 데 사용됩니다.
Vincenty 공식
- 정의: Vincenty 공식은 1975년 Thaddeus Vincenty에 의해 개발된 방법으로, 지구상의 두 점 사이의 최단 거리를 계산합니다. 이 공식은 WGS-84 같은 타원체 모델을 사용하여 높은 정확도로 거리를 측정합니다.
- 작동 원리: Vincenty 공식은 반복적인 방법을 사용하여 두 지점 사이의 거리를 계산합니다. 이 과정에서, 공식은 타원체의 장반경, 편평률 및 두 지점의 위도와 경도를 고려합니다.
- 정확성과 한계: Vincenty 공식은 대부분의 경우 매우 정확한 결과를 제공하지만, 매우 먼 거리(특히 거의 안티포달 지점, 즉 서로 반대편에 위치한 지점)에서는 수렴하지 않을 수 있습니다. 이러한 경우에는 다른 방법(예: Haversine 공식)을 사용하는 것이 좋습니다.
Vincenty 공식은 지리학, 항해, GPS 시스템 설계 및 위치 기반 서비스 개발과 같은 다양한 분야에서 중요한 도구로 사용됩니다. 이 공식은 높은 정확도를 요구하는 응용 프로그램에 특히 유용합니다.
package main
import (
"math"
)
// WGS-84 ellipsoid parameters
const (
a = 6378137.0 // semi-major axis
f = 1 / 298.257223563 // flattening
b = (1 - f) * a // semi-minor axis
)
// toRadians converts degrees to radians
func toRadians(deg float64) float64 {
return deg * math.Pi / 180
}
// vincentyDistance calculates the distance between two points on the Earth
func vincentyDistance(lat1, lon1, lat2, lon2 float64) float64 {
// Convert latitude and longitude from degrees to radians
U1 := math.Atan((1 - f) * math.Tan(toRadians(lat1)))
U2 := math.Atan((1 - f) * math.Tan(toRadians(lat2)))
L := toRadians(lon2 - lon1)
Lambda := L
var sinSigma, cosSigma, sigma, sinAlpha, cos2Alpha, cos2SigmaM, C float64
for i := 0; i < 1000; i++ {
sinLambda := math.Sin(Lambda)
cosLambda := math.Cos(Lambda)
sinSigma = math.Sqrt(math.Pow(math.Cos(U2)*sinLambda, 2) +
math.Pow(math.Cos(U1)*math.Sin(U2)-math.Sin(U1)*math.Cos(U2)*cosLambda, 2))
cosSigma = math.Sin(U1)*math.Sin(U2) + math.Cos(U1)*math.Cos(U2)*cosLambda
sigma = math.Atan2(sinSigma, cosSigma)
sinAlpha = math.Cos(U1) * math.Cos(U2) * sinLambda / sinSigma
cos2Alpha = 1 - sinAlpha*sinAlpha
cos2SigmaM = cosSigma - 2*math.Sin(U1)*math.Sin(U2)/cos2Alpha
C = f / 16 * cos2Alpha * (4 + f*(4-3*cos2Alpha))
LambdaPrev := Lambda
Lambda = L + (1-C)*f*sinAlpha*
(sigma+C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)))
// Check for convergence
if math.Abs(Lambda-LambdaPrev) < 1e-12 {
break
}
}
uSquared := cos2Alpha * (a*a - b*b) / (b * b)
A := 1 + uSquared/16384*(4096+uSquared*(-768+uSquared*(320-175*uSquared)))
B := uSquared / 1024 * (256 + uSquared*(-128+uSquared*(74-47*uSquared)))
deltaSigma := B * sinSigma * (cos2SigmaM + B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)))
// Calculate the distance
s := b * A * (sigma - deltaSigma)
return s // distance in meters
}
func main() {
lat1, lon1 := 37.7749, -122.4194 // San Francisco
lat2, lon2 := 34.0522, -118.2437 // Los Angeles
distance := vincentyDistance(lat1, lon1, lat2, lon2)
println("Distance:", distance, "meters")
}
위 소스 코드를 살펴보면 vincentyDistance함수에 두 지점의 위경도 값을 실수 형태로 전달하면, 두 지점 사이의 거리를 미터(m) 단위로 반환합니다.
따라서 위 소스코드를 실행하면 샌프란시스코와 로스앤젤레스 사이 거리를 +5.590423e+005 meters로 반환합니다.
(본 글은 지금은 운영하지 않는 과거 저의 블로그에 있었던 글입니다.)
