将纬度和经度转换为 3D 空间中的点

Posted 2023-03-05

【中文标题】将纬度和经度转换为 3D 空间中的点

【英文标题】：Convert Latitude and Longitude to point in 3D space

【发布时间】：2012-05-06 20:07:34

【问题描述】：

我需要将纬度和经度值转换为 3 维空间中的一个点。我已经尝试了大约 2 个小时，但我没有得到正确的结果。

Equirectangular 坐标来自openflights.org。我尝试了几种 cos 和 sin 的组合，但结果看起来并不像我们心爱的小地球。

---

在下面，您可以看到应用Wikipedia 建议的转换的结果。我想可以从上下文中猜出c4d.Vector 是什么。

def llarToWorld(latit, longit, altid, rad):
    x = math.sin(longit) * math.cos(latit)
    z = math.sin(longit) * math.sin(latit)
    y = math.cos(longit)
    v = c4d.Vector(x, y, z)
    v = v * altid + v * rad
    return v

红色：X，绿色：Y，蓝色：Z

确实可以识别北美和南美，尤其是墨西哥湾周围的土地。但是，它看起来有点压扁，有点放错地方了..

---

由于结果看起来有些旋转，我想，我尝试交换纬度和经度。但是这个结果有点尴尬。

def llarToWorld(latit, longit, altid, rad):
    temp = latit
    latit = longit
    longit = temp
    x = math.sin(longit) * math.cos(latit)
    z = math.sin(longit) * math.sin(latit)
    y = math.cos(longit)
    v = c4d.Vector(x, y, z)
    v = v * altid + v * rad
    return v

---

这是没有转换值的结果。

def llarToWorld(latit, longit, altid, rad):
    return c4d.Vector(math.degrees(latit), math.degrees(longit), altid)

---

问题：如何正确转换经纬度？

---

解决方案

感谢 TreyA，我在 mathworks.com 上找到了 this 页面。执行它的代码如下：

def llarToWorld(lat, lon, alt, rad):
    # see: http://www.mathworks.de/help/toolbox/aeroblks/llatoecefposition.html
    f  = 0                              # flattening
    ls = atan((1 - f)**2 * tan(lat))    # lambda

    x = rad * cos(ls) * cos(lon) + alt * cos(lat) * cos(lon)
    y = rad * cos(ls) * sin(lon) + alt * cos(lat) * sin(lon)
    z = rad * sin(ls) + alt * sin(lat)

    return c4d.Vector(x, y, z)

实际上，我切换了y 和z，因为那时地球是旋转的，但是，它起作用了！结果是这样的：

【问题讨论】：

altid 是海拔高度，但rad 是什么？那是地球的半径吗？ altid 和 rad 是否使用相同的单位（英尺）？如果您只使用半径（即仅v = v * rad）怎么办？

还可以查看谷歌搜索“lla to ecef” - 以地球为中心的纬度/经度/高度，固定地球。

@TreyA 完美，谢谢！发现这个：mathworks.de/help/toolbox/aeroblks/llatoecefposition.html 这是正确的公式。 :) 如果您想要代表，您可以将您的评论作为答案，这样我也可以将我的问题标记为已回答。

【参考方案1】：

我已经重新格式化了之前在这里提到的代码，但更重要的是，您遗漏了 Niklas R 提供的链接中提到的一些方程式

def LLHtoECEF(lat, lon, alt):
    # see http://www.mathworks.de/help/toolbox/aeroblks/llatoecefposition.html

    rad = np.float64(6378137.0)        # Radius of the Earth (in meters)
    f = np.float64(1.0/298.257223563)  # Flattening factor WGS84 Model
    cosLat = np.cos(lat)
    sinLat = np.sin(lat)
    FF     = (1.0-f)**2
    C      = 1/np.sqrt(cosLat**2 + FF * sinLat**2)
    S      = C * FF

    x = (rad * C + alt)*cosLat * np.cos(lon)
    y = (rad * C + alt)*cosLat * np.sin(lon)
    z = (rad * S + alt)*sinLat

    return (x, y, z)

比较输出：查找加利福尼亚州洛杉矶的 ECEF（34.0522，-118.40806，0 海拔） 我的代码： X = -2516715.36114 米或 -2516.715 公里 Y = -4653003.08089 米或 -4653.003 公里 Z = 3551245.35929 米或 3551.245 公里

您的代码： X = -2514072.72181 米或-2514.072 公里 Y = -4648117.26458 米或-4648.117 公里 Z = 3571424.90261 米或 3571.424 公里

虽然在您的地球自转环境中，您的函数会生成正确的地理区域以供显示，但它不会给出正确的 ECEF 等效坐标。正如您所看到的，一些参数的差异高达 20 KM，这是一个相当大的错误。

扁平化因子，f 取决于您为转换假设的模型。典型，型号为WGS 84；但是，还有其他型号。

就我个人而言，我喜欢在海军研究生院使用this link 来检查我的转化情况。

【讨论】：

这可行，但您首先需要将角度转换为弧度。纬度/经度值通常以度数形式提供。

正如@TomislavMuic 所指出的，您必须在使用上述 sn-p 之前将 degree lat lon 转换为弧度。为此，只需使用 math.radians() 函数。

【参考方案2】：

你没有按照***的建议去做。请仔细阅读。

他们说：

x = r cos(phi) sin(theta)
y = r sin(phi) sin(theta)
z = r cos(theta)

然后：

theta == latitude
phi == longitude

在你的情况下，r = 半径 + 高度

所以你应该使用：

r = radius + altitude
x = r cos(long) sin(lat)
y = r sin(long) sin(lat)
z = r cos(lat)

请注意，最后一个条目是cos(lat)（您使用的是经度）。

【讨论】：

不幸的是，这与像素 #2 显示的结果相同.. :(

不，不是。 lat/long 仅交换 z。你把它们都换了。

嗯，不完全一样，但和第二张图差不多。交换纬度和经度时，它与第一张图片非常相似。我通过搜索lla to ecef（感谢 TreyA）找到的链接非常完美。见：mathworks.de/help/toolbox/aeroblks/llatoecefposition.html

应该是x = r cos(long) cos(lat) y = r sin(long) cos(lat) z = r sin(lat)。请注意 lat 与 xy 平面的角度，而不是与***建议的 z 轴的角度。

【参考方案3】：

正如 TreyA 所说，LLA to ECEF 是解决方案。见http://www.mathworks.de/help/toolbox/aeroblks/llatoecefposition.html