将纬度和经度转换为三维空间中的点

25

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

Equirectangular坐标来自openflights.org。我尝试了几种cossin的组合,但结果从未像我们心爱的地球那样。


在下面,您可以看到应用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

enter image description here

红色: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)

enter image description here


Solution

Thanks to TreyA, I found this page on mathworks.com, which provides a code snippet for converting longitude and latitude. Here's the code:

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)

实际上,我交换了yz因为那时地球被旋转了,不过它可以工作!这就是结果:

输入图像描述


altid 是海拔高度,但 rad 是什么?是地球半径吗?altidrad 是否使用相同的单位(英尺)?如果只使用半径会怎样(即只有 v = v * rad)? - Omri Barel
还可以搜索“lla to ecef” - 将纬度/经度/高度转换为地心坐标系。 - TreyA
1
@TreyA 太好了,谢谢!找到了这个链接:http://www.mathworks.de/help/toolbox/aeroblks/llatoecefposition.html 这就是正确的公式。 :) 如果你想要声望分,你可以把你的评论变成一个答案,这样我也可以标记我的问题已经得到解答了。 - Niklas R
3个回答

18

我重新格式化了之前提到的代码,但更重要的是你忽略了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)

比较输出:查找洛杉矶,加利福尼亚州(34.0522,-118.40806,0海拔)的ECEF坐标
我的代码:
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公里,相当大的误差。

扁率因子f取决于您假设的转换模型。 典型模型是WGS 84; 但是,还有其他模型。

就我个人而言,我喜欢使用这个链接来对我的转换进行健全性检查。


这个可以运行,但是你需要先将角度转换为弧度。而经纬度值通常以度数提供。 - Tomislav Muic
正如@TomislavMuic所指出的那样,在使用上述代码片段之前,您必须将degree纬度和经度转换为弧度。要做到这一点,只需使用math.radians()函数即可。 - Harshdeep

15

您没有按照维基百科的建议去做。请再仔细阅读一遍。

他们说:

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显示的结果相同.. :( - Niklas R
不是这样的。只有 z 的经纬度被交换了。你把它们全部交换了。 - andrew cooke
嗯,虽然不完全一样,但与第二张图片非常相似。当经纬度互换时,与第一张图片非常相近。我通过搜索“lla to ecef”找到的链接(感谢TreyA)非常完美。请查看:http://www.mathworks.de/help/toolbox/aeroblks/llatoecefposition.html - Niklas R
3
应该这样写:x = r cos(long) cos(lat) y = r sin(long) cos(lat) z = r sin(lat)。注意,lat是相对于 xy 平面的角度,而不是维基百科所建议的相对于 z 轴的角度。 - sachinruk

4

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接