개인적으로 부동산 분석 토이 프로젝트를 진행하다가 두 GPS 좌표간 거리를 구해야 할 일이 생겼습니다. 해당 경험이 전무한 저로서는 (무식하면 용감하다고) 구글링을 통해 찾게 된 코드를 바로 써먹을 생각에 신이 나 있었는데요, 그 코드는 아래와 같습니다.

from math import sin, cos, sqrt, atan2, radians

# approximate radius of earth in km
R = 6373.0

lat1 = radians(52.2296756)
lon1 = radians(21.0122287)
lat2 = radians(52.406374)
lon2 = radians(16.9251681)

dlon = lon2 - lon1
dlat = lat2 - lat1

a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
c = 2 * atan2(sqrt(a), sqrt(1 - a))

distance = R * c

print("Result:", distance)
print("Should be:", 278.546, "km")

(한가지 참고 목적으로 짚고 넘어갈 점은, 파이썬 math 모듈의 삼각함수는 degree가 아닌 radian을 사용한다는 점입니다.)

그런데 여기서 Kurt Peek의 날카로운 답변을 보게 됩니다. 지구는 완전한 구가 아닌 회전타원체에 가깝고 이로 인해 0.5%에 상당하는 오차가 생길 수 있다는 점 때문에 Haversine distance 대신 WGS-84 좌표 시스템을 사용하는 Vincenty distance를 사용할 것을 추천했습니다.

The answers above are based on the Haversine formula, which assumes the earth is a sphere, which results in errors of up to about 0.5% (according to help(geopy.distance)). Vincenty distance uses more accurate ellipsoidal models such as WGS-84, and is implemented in geopy.

다만 최근의 geopy documentation을 보면, Vincenty(1975) 역시 deprecate 되었고, 현재는 Karney(2013)를 디폴트 알고리즘으로 하고 있다는 점을 알게 되었습니다.

The geodesic distance is the shortest distance on the surface of an ellipsoidal model of the earth. The default algorithm uses the method is given by Karney (2013) (geodesic); this is accurate to round-off and always converges. An older deprecated method due to Vincenty (1975) (vincenty) is also available; this is only accurate to 0.2 mm and the distance calculation fails to converge for nearly antipodal points.

from geopy import distance
newport_ri = (41.49008, -71.312796)
cleveland_oh = (41.499498, -81.695391)
print(distance.distance(newport_ri, cleveland_oh).miles)

wellington = (-41.32, 174.81)
salamanca = (40.96, -5.50)
print(distance.distance(wellington, salamanca).km)