我有一些以地球为中心的坐标点作为纬度和经度 (WGS-84)。
如何将它们转换为原点位于地球中心的笛卡尔坐标 (x,y,z)？
【问题讨论】：
标签： mapping geometry geospatial
这是我找到的答案：
只是为了使定义完整，在笛卡尔坐标系中：
转换为：
x = R * cos(lat) * cos(lon)
y = R * cos(lat) * sin(lon)
z = R *sin(lat)
其中 R 是 the approximate radius of earth（例如 6371 公里）。
如果您的三角函数需要弧度（他们可能会这样做），您需要先将经度和纬度转换为弧度。您显然需要十进制表示，而不是度\分\秒（请参阅e.g. here 关于转换）。
反向转换公式：
lat = asin(z / R)
lon = atan2(y, x)
asin 当然是反正弦。 read about atan2 in wikipedia。不要忘记将弧度转换回度数。
This page为此提供了c#代码（注意它与公式有很大不同），还有一些解释和漂亮的图表说明为什么这是正确的，
【讨论】：
我最近使用 WGS-84 数据上的“Haversine 公式”，它是“Haversines 定律”的衍生物，结果非常令人满意。
是的，WGS-84 假设地球是一个椭球体，但我相信使用“Haversine 公式”之类的方法只会得到大约 0.5% 的平均误差，这在你的情况下可能是可接受的误差量。除非您谈论的是几英尺的距离，否则您总会有一些错误，即使在理论上地球也存在曲率......如果您需要更严格的 WGS-84 兼容方法，请查看“Vincenty 公式。 "
我了解 starblue 的来源，但好的软件工程通常需要权衡取舍，因此这完全取决于您所做工作所需的准确性。例如，从“曼哈顿距离公式”计算的结果与“距离公式”的结果相比，在某些情况下可能会更好，因为它的计算成本更低。想想“哪个点最近？”不需要精确距离测量的场景。
关于“Haversine 公式”，它很容易实现并且很好，因为它使用“球面三角”而不是基于“余弦定律”的基于二维三角函数的方法，因此你会得到一个很好的准确性与复杂性之间的平衡。
一位名叫Chris Veness的绅士有一个很棒的website，它解释了您感兴趣的一些概念并演示了各种程序化实现；这也应该回答您的 x/y 转换问题。
【讨论】：
在 python3.x 中可以使用：
# Converting lat/long to cartesian
import numpy as np
def get_cartesian(lat=None,lon=None):
lat, lon = np.deg2rad(lat), np.deg2rad(lon)
R = 6371 # radius of the earth
x = R * np.cos(lat) * np.cos(lon)
y = R * np.cos(lat) * np.sin(lon)
z = R *np.sin(lat)
return x,y,z
【讨论】：
将GPS（WGS84) 转换为笛卡尔坐标的理论
https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates
以下是我正在使用的：
我附上了一个我写的VB代码：
Imports System.Math
'Input GPSLatitude is WGS84 Latitude,h is altitude above the WGS 84 ellipsoid
Public Function GetSphericalLatitude(ByVal GPSLatitude As Double, ByVal h As Double) As Double
Dim A As Double = 6378137 'semi-major axis
Dim f As Double = 1 / 298.257223563 '1/f Reciprocal of flattening
Dim e2 As Double = f * (2 - f)
Dim Rc As Double = A / (Sqrt(1 - e2 * (Sin(GPSLatitude * PI / 180) ^ 2)))
Dim p As Double = (Rc + h) * Cos(GPSLatitude * PI / 180)
Dim z As Double = (Rc * (1 - e2) + h) * Sin(GPSLatitude * PI / 180)
Dim r As Double = Sqrt(p ^ 2 + z ^ 2)
Dim SphericalLatitude As Double = Asin(z / r) * 180 / PI
Return SphericalLatitude
End Function
请注意h 的高度高于WGS 84 ellipsoid。
通常GPS 会给我们H 高于MSL 的高度。
必须使用 geopotential 模型 EGM96 将 MSL 高度转换为 h 高于 WGS 84 ellipsoid 的高度（Lemoine et al, 1998）。
这是通过对空间分辨率为 15 角分的大地水准面高度文件的网格进行插值来完成的。
或者如果你有一定的水平专业GPS有海拔H（msl，高度高于平均海平面）和UNDULATION之间的关系geoid 和 ellipsoid (m) 从内部表中选择的基准输出。你可以得到h = H(msl) + undulation
通过笛卡尔坐标到 XYZ：
x = R * cos(lat) * cos(lon)
y = R * cos(lat) * sin(lon)
z = R *sin(lat)
【讨论】：
如果您关心基于椭球而不是球体获取坐标，请查看Geographic_coordinate_conversion - 它提供了公式。 GEodetic Datum 具有转换所需的 WGS84 常量。
那里的公式还考虑了相对于参考椭球表面的高度（如果您从 GPS 设备获取高度数据，则很有用）。
【讨论】：
为什么要实施已经实施和测试证明的东西？
例如，C# 具有 NetTopologySuite，它是 JTS 拓扑套件的 .NET 端口。
具体来说，您的计算存在严重缺陷。地球不是一个完美的球体，地球半径的近似值可能无法用于精确测量。
如果在某些情况下使用自制函数是可以接受的，那么 GIS 就是一个很好的例子，在该领域中更倾向于使用可靠的、经过测试证明的库。
【讨论】：
Coordinate[] coordinates = new Coordinate[3];
coordinates[0] = new Coordinate(102, 26);
coordinates[1] = new Coordinate(103, 25.12);
coordinates[2] = new Coordinate(104, 16.11);
CoordinateSequence coordinateSequence = new CoordinateArraySequence(coordinates);
Geometry geo = new LineString(coordinateSequence, geometryFactory);
CoordinateReferenceSystem wgs84 = DefaultGeographicCRS.WGS84;
CoordinateReferenceSystem cartesinaCrs = DefaultGeocentricCRS.CARTESIAN;
MathTransform mathTransform = CRS.findMathTransform(wgs84, cartesinaCrs, true);
Geometry geo1 = JTS.transform(geo, mathTransform);
【讨论】：
java.lang.IllegalArgumentException: dimension must be <= 3
我在 Python 中创建了一个函数，考虑到地球不是一个完美的球体这一事实。参考文献在 cmets 中：
# this function converts latitude,longitude and height above sea level
# to earthcentered xyx coordinates in wgs84, lat and lon in decimal degrees
# e.g. 52.724156(West and South are negative), heigth in meters
# for algoritm see https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates
# for values of a and b see https://en.wikipedia.org/wiki/Earth_radius#Radius_of_curvature
from math import *
def latlonhtoxyzwgs84(lat,lon,h):
a=6378137.0 #radius a of earth in meters cfr WGS84
b=6356752.3 #radius b of earth in meters cfr WGS84
e2=1-(b**2/a**2)
latr=lat/90*0.5*pi #latitude in radians
lonr=lon/180*pi #longituede in radians
Nphi=a/sqrt(1-e2*sin(latr)**2)
x=(Nphi+h)*cos(latr)*cos(lonr)
y=(Nphi+h)*cos(latr)*sin(lonr)
z=(b**2/a**2*Nphi+h)*sin(latr)
return([x,y,z])
【讨论】：
你可以在 Java 上这样做。
public List<Double> convertGpsToECEF(double lat, double longi, float alt) {
double a=6378.1;
double b=6356.8;
double N;
double e= 1-(Math.pow(b, 2)/Math.pow(a, 2));
N= a/(Math.sqrt(1.0-(e*Math.pow(Math.sin(Math.toRadians(lat)), 2))));
double cosLatRad=Math.cos(Math.toRadians(lat));
double cosLongiRad=Math.cos(Math.toRadians(longi));
double sinLatRad=Math.sin(Math.toRadians(lat));
double sinLongiRad=Math.sin(Math.toRadians(longi));
double x =(N+0.001*alt)*cosLatRad*cosLongiRad;
double y =(N+0.001*alt)*cosLatRad*sinLongiRad;
double z =((Math.pow(b, 2)/Math.pow(a, 2))*N+0.001*alt)*sinLatRad;
List<Double> ecef= new ArrayList<>();
ecef.add(x);
ecef.add(y);
ecef.add(z);
return ecef;
}
【讨论】：