구에서 좌표사이의 거리 구하기

직교좌표계에서 내적을 구하는 건 매우 간단하다. 내적은 단위원 위의 두 벡터에 대해서는 두 벡터 사이의 사이각의 코사인. 최단거리이기 때문에 180도 이상인지 이하인지는 중요하지 않다. 그래서 내적의 역코사인에 구의 반지름을 곱하면 답.

 

void GetCartesian(double Lat, double Long, double *x, double *y, double *z)
{
/*
    유사구면좌표계를 직교좌표계로 변환.
    http://ko.wikipedia.org/wiki/%EA%B5%AC%EB%A9%B4_%EC%A2%8C%ED%91%9C%EA%B3%84
 */
    *x = cos(Lat)*cos(Long);
    *y = cos(Lat)*sin(Long);
    *z = sin(Lat);
}
double Calc1(double Lat1, double Long1, double Lat2, double Long2)
{
/*
    구면위의 두 점 사이의 사이각(theta<=PI)만 필요하다. 두 좌표의 내적으로부터 구한다.
*/
    double dDistance;
    const double kEarthRadiusKms = 6376.5;
    double Lat1InRad  = Lat1 * (PI / 180.0);
    double Long1InRad = Long1 * (PI / 180.0);
    double Lat2InRad  = Lat2 * (PI / 180.0);
    double Long2InRad = Long2 * (PI / 180.0);

    double x1, y1, z1;
    double x2, y2, z2;
    double theta;

    GetCartesian(Lat1InRad, Long1InRad, &x1, &y1, &z1);
    GetCartesian(Lat2InRad, Long2InRad, &x2, &y2, &z2);

    theta = acos(x1*x2 + y1*y2 + z1*z2); // 내적 = cos(theta)

    // Distance.
    dDistance = kEarthRadiusKms * theta;

    return dDistance;
}

from typing import Tuple
from math import cos, sin, PI, acos


EARTH_RADIUS_KM = 6376.5


def to_cartesian(lat: float, lon: float) -> Tuple[float, float, float]:
    x = cos(lat) * cos(lon)
    y = cos(lat) * sin(lon)
    z = sin(lat)
    return x, y, z


def degree_to_radian(deg: float) -> float:
    return deg * (PI / 180)


def distance_on_sphere(pos1: Tuple[float, float], pos2: Tuple[float, float]) -> float:

    pos1_radian = list(map(degree_to_radian, pos1))
    pos2_radian = list(map(degree_to_radian, pos2))

    pos1_xyz = to_cartesian(*pos1_radian)
    pos2_xyz = to_cartesian(*pos2_radian)

    theta = acos(sum(e1 * e2 for e1, e2 in zip(pos1_xyz, pos2_xyz)))

    distance = EARTH_RADIUS_KM * theta

    return distance