如何使用Python比较两个球面点的大圆距离和欧几里得距离？

Question

我正在尝试检查当您使用欧几里得距离计算地球上两点之间的距离而不是使用大圆距离（gcd）时引入的错误。 我有两个由其纬度和经度定义的点。 我使用了Python geopy框架来计算大圆距离。这是计算gcd的代码：

def measure(self, a, b):
        a, b = Point(a), Point(b)

        lat1, lng1 = radians(degrees=a.latitude), radians(degrees=a.longitude)
        lat2, lng2 = radians(degrees=b.latitude), radians(degrees=b.longitude)

        sin_lat1, cos_lat1 = sin(lat1), cos(lat1)
        sin_lat2, cos_lat2 = sin(lat2), cos(lat2)

        delta_lng = lng2 - lng1
        cos_delta_lng, sin_delta_lng = cos(delta_lng), sin(delta_lng)

        d = atan2(sqrt((cos_lat2 * sin_delta_lng) ** 2 +
                       (cos_lat1 * sin_lat2 -
                        sin_lat1 * cos_lat2 * cos_delta_lng) ** 2),
                  sin_lat1 * sin_lat2 + cos_lat1 * cos_lat2 * cos_delta_lng)

        return self.RADIUS * d

有两个要点：

p1=[39.8616,-75.0748], p2=[-7.30933,112.76]

其中

gcd = 78.8433004543197 klm

使用geopy库中的great_circle(p1,p2).kilometers函数计算出了这两点之间的大圆距离（gcd）。

接着，我使用以下公式将这两个点转换成笛卡尔坐标：

  def spherical_to_cartesian(r,la,lo):
       x=r*np.sin(90-la)*np.cos(lo)
       y=r*np.sin(90-la)*np.sin(lo)
       z=r*np.cos(90-la)
       return (x,y,z)

其中r=6372.795，导致以下笛卡尔坐标：

p1=[ -765.81579368,  -256.69640558,  6321.40405587], 
p2=[480.8302149,-168.64726394,-6352.39140142]

然后通过输入: np.linalg.norm(p2-p1)，我得到了它们的欧几里得范数为1103.4963114787836，与从GCD中得到的约78公里相比似乎不太合理。我是否推断错了什么？

- curious

---

你确定你在正确的单位进行比较吗？ - jonrsharpe

不确定这是否是您唯一的错误，但 np.sin 和 np.cos 是以弧度为单位运算的... - John Greenall

GDC看起来非常小。你有两个点相隔很多度。我认为位置上的1度差异将是100公里。 - Salix alba

你是对的！问题在于我计算笛卡尔坐标的方式。我应该先转换为弧度。现在我得到了78.842797622886096，非常接近最大公约数。 - curious

如果你想知道的话，从新泽西州收费站到印度尼西亚苏拉巴亚的某个地方... - Nick T

Answer 1

Python在math模块中提供了两个函数；radians将角度转换为弧度，degrees将弧度转换为角度。

方法sin()返回x的正弦值，单位为弧度。

import math
def spherical_to_cartesian(r,la,lo):
   rlo = math.radians(lo)
   rla = math.radians(90-la)
   x=r*np.sin(rla)*np.cos(rlo)
   y=r*np.sin(rla)*np.sin(rlo)
   z=r*np.cos(rla)
   return (x,y,z)

- sabbahillel