Skyfield中的地球形状似乎有问题 - 我的Python代码正确吗？

Question

将地球中心到各种 (纬度，经度) 位置的距离映射使用Skyfield，显示与纬度有关但与经度无关（亚毫米级）。 这可能是包中记录的近似值，我的脚本中的错误或其他原因。 我在这里做错了什么吗？（除了当然使用喷气式飞机）

import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import load, now

data  = load('de421.bsp')
earth = data['earth']
jd    = now()

epos = earth.at(jd).position.km

lats = np.linspace( -90,  90, 19)
lons = np.linspace(-180, 180, 37)
LATS, LONS = np.meshgrid(lats, lons)

s = LATS.shape

points = zip(LATS.flatten(), LONS.flatten())

rr = []
for point in points:
    la, lo = point
    pos = earth.topos(la, lo).at(jd).position.km
    r   = np.sqrt( ((pos-epos)**2).sum() )
    rr.append(r)

surf = np.array(rr).reshape(s)

extent = [lons.min(), lons.max(), lats.min(), lats.max()]

plt.figure()
plt.imshow(surf.T, origin='lower', extent=extent)
plt.colorbar()
plt.title('uhoh topo')
plt.savefig('uhoh topo')
plt.show()

作为交叉检验，我随机尝试了一些纬度相同的位置：

pe = earth.at(jd).position.km
for i in range(10):

    lon1, lon2 = 360.*np.random.random(2)-180
    lat = float(180.*np.random.random(1)-90.)

    p1  = earth.topos(lat, lon1).at(jd).position.km
    p2  = earth.topos(lat, lon2).at(jd).position.km

    r1   = np.sqrt( ((p1-pe)**2).sum() )
    r2   = np.sqrt( ((p2-pe)**2).sum() )

    print lat, lon1, lon2, r2-r1

并得到以下结果（第四列显示微米的差异）：

45.8481950437 55.9538249618 115.148786114 1.59288902069e-08
-72.0821405192 4.81264755835 172.783338907 2.17096385313e-09
51.6126938075 -54.5670258363 -134.888403816 2.42653186433e-09
2.92691713179 -178.553103457 134.648099589 1.5916157281e-10
-78.7376163827 -55.0684703115 125.714124504 -6.13908923697e-10
48.5852207923 -169.061708765 35.5374862329 7.60337570682e-10
42.3767785876 130.850223447 -111.520896867 -1.62599462783e-08
11.2951212126 -60.0296460731 32.8775784623 6.91579771228e-09
18.9588262131 71.3414406837 127.516370219 -4.84760676045e-09
-31.5768658495 173.741960359 90.3715297869 -6.78483047523e-10

- uhoh

---

实际问题是什么？从大约6357公里到6378公里的纬度变化看起来合理（例如，请参见https://en.wikipedia.org/wiki/Figure_of_the_Earth），因此您在输出中质疑哪些数字？ - Warren Weckesser

距离应该随着经度变化而变化 - 肯定比几十微米更多。扁球体只是对真实形状的一阶近似。这甚至没有考虑高程或大地水准面的形状。我试图理解这是否是软件包中的隐式近似，还是一个错误，或者我没有正确使用软件包。 - uhoh

经度变化虽然很小但不为零，可能是由于四舍五入误差（这些是天文距离），因此有可能使用了椭球近似。但我不知道是否是这种情况。人们进行GPS测量和卫星定位以毫米分辨率（有时还要准确性！）-如果是这种情况，那将是一个相当大的近似。 - uhoh

Answer 1

< p > topos() 方法包括一个 elevation_m=0.0 ，您没有指定。因此，在每次调用中，它都采用默认值为海平面以上0.0米。而“海平面”的定义不随经度变化 - 在给定纬度上，海平面与地球中心的距离始终相同。

我不确定为什么您将微米级误差称为“相当大的近似值”，但当您减去大约为1 au的两个距离时，您确实遇到了机器64位浮点数学的有限精度 - 您看到了最后一位或两位的舍入误差。

- Brandon Rhodes

---

我检查了earth.topos.__doc__ - 应该输入help(earth.topos)，其中包含高程选项。在线文档显示lat和lon参数，但没有高程m。我知道文档正在进行中 - 也许只需将高程选项添加到卫星示例中，并使用一个不在海平面上（丹佛而不是波士顿？）的城市。此外，我认为即使elevation = 0海平面表面也会是一个大地水准面 - 重力等势面，并且具有2D结构。现在我看到WGS84双轴椭球体对于大地测量目的来说是“海平面”。谢谢 - Skyfield非常棒！ - uhoh

“相当大的近似值”是指仅使用椭球体来计算近地和天文事件，这是对@WarrenWeckesser的评论做出的回应，在我还不知道有高程选项可用之前。你肯定会同意，不包括高程将是一个很大的近似值！顺便说一句，如果你能看一下这个问题，那就太好了。再次感谢！ - uhoh