Converting from longitude\latitude to Cartesian coordinates

我有一些以地球为中心的坐标点作为经纬度(WGS-84)。

如何将它们转换为笛卡尔坐标系(x，y，z) ，原点在地球的中心？

331233 次浏览
小开

这是我找到的答案:

只是为了让定义更完整，在笛卡儿坐标系:

  • X 轴经过长纬度(0,0) ，所以经度0与赤道相交;
  • the y-axis goes through (0,90);
  • and the z-axis goes through the poles.

转换是:

x = R * cos(lat) * cos(lon)


y = R * cos(lat) * sin(lon)


z = R *sin(lat)

其中 R 为 地球的近似半径(例如6371公里)。

如果你的三角函数期望得到弧度(他们可能会这么做) ，你需要先把你的经纬度转换成弧度。你显然需要一个小数，而不是度分钟(见 例如这里关于转换)。

返回转换公式:

   lat = asin(z / R)
lon = atan2(y, x)

当然是弧正弦。 read about atan2 in wikipedia。别忘了把弧度转换回度数。

这个页面 给出了这个问题的 c # 代码(请注意，它与公式非常不同) ，还有一些解释和为什么这是正确的图表,

小开

为什么要实现已经实现并经过测试验证的东西？

例如，C # 有 NetTopologySuite，它是 JTS 拓扑套件的.NET 端口。

具体来说，你的计算有一个严重的缺陷。地球不是一个完美的球体，地球半径的 approximation可能无法精确测量。

如果在某些情况下，使用自制函数是可以接受的，那么 GIS 就是一个很好的例子，在这个领域中，人们更倾向于使用可靠的、经过测试验证的库。

小开
最佳答案

我最近使用 “半正矢公式”在 WGS-84数据上，这是“哈佛定律”的衍生物，结果非常令人满意。

是的，WGS-84假设地球是一个椭球体，但是我相信使用“半正矢公式”这样的方法，你只能得到0.5% 的平均误差，在你的情况下，这可能是一个可以接受的误差量。你总是会有一些误差，除非你说的是几英尺的距离，即使那样，理论上地球是弯曲的... 如果你需要一个更严格的 WGS-84兼容的方法检查“文森蒂公式”

我知道 星蓝是从哪里来的，但是好的软件工程通常是关于权衡的，所以这完全取决于您所做的事情所需要的准确性。例如，从“曼哈顿距离公式”计算的结果与从“距离公式”计算的结果在某些情况下可以更好，因为它的计算成本较低。想想“哪个点最近?”不需要精确测量距离的场景。

关于“半正矢公式”，它很容易实现，也很好，因为它使用的是“球面三角学”，而不是基于二维三角函数的“余弦定律”方法，因此你可以在精度和复杂性之间取得一个很好的平衡。

一位名为 克里斯 · 维内斯的绅士有一个很棒的 网站，它解释了您感兴趣的一些概念，并演示了各种编程实现; 这应该也可以回答您的 x/y 转换问题。

小开

如果你关心基于椭球体而不是球体的坐标，看看 地理坐标转换-它给出了公式。大地测量基准具有转换所需的 WGS84常量。

这里的公式还考虑了相对于参考椭球体表面的高度(如果你正在从 GPS 设备获取高度数据，这个公式很有用)。

小开 
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);
小开

将 GPS（WGS84)转换为 < strong > < em > 笛卡尔坐标系的理论 Https://en.wikipedia.org/wiki/geographic_coordinate_conversion#from_geodetic_to_ecef_coordinates

以下是我正在使用的:

  • 经度在 GPS (WGS84)和笛卡尔坐标是相同的。
  • 需要用 WGS 84椭球参数转换的纬度半长轴为6378137米
  • 压扁的倒数是298.257223563。

我附上了一份 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高度。 通过使用 地理位置型 EGM96(Lemoine et al, 1998) ，MSL高度必须转换为高于 WGS 84 ellipsoid的 h高度。
这是通过以15弧分的空间分辨率插值大地水准面高度文件的网格来完成的。

或者如果你有一些级别的 很专业，GPS有高度 H(高于平均海平面)和 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)
小开

项目4软件提供了一个命令行程序，可以进行转换，例如。

LAT=40
LON=-110
echo $LON $LAT | cs2cs +proj=latlong +datum=WGS84 +to +proj=geocent +datum=WGS84

它还提供了一个 C API。特别是，函数 pj_geodetic_to_geocentric将不必首先设置投影对象就可以进行转换。

小开

你可以在 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;




}
小开

在 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
小开

我用 Python 编写了一个函数，它考虑到地球不是一个完美的球体这一事实。参考资料见评论:

    # 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])

提交答案