114

我有一些以地球为中心的坐标点，以纬度和经度（WGS-84）形式给出。

如何将它们转换为原点位于地球中心的笛卡尔坐标（x，y，z）？

4

10 回答 10

153

这是我找到的答案：

只是为了使定义完整，在笛卡尔坐标系中：

  • x轴经过long,lat(0,0)，所以经度0与赤道相交；
  • y 轴穿过 (0,90)；
  • z轴穿过两极。

转换是：

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)

asin 当然是反正弦。在维基百科中阅读 atan2。不要忘记将弧度转换回度数。

此页面为此提供了 c# 代码（请注意，它与公式有很大不同），还有一些解释和很好的图表说明为什么这是正确的，

于 2009-07-26T20:02:13.460 回答
53

我最近在 WGS-84 数据上使用“Haversine 公式”做了类似的事情，这是“Haversines 定律”的衍生物，结果非常令人满意。

是的，WGS-84 假设地球是一个椭球体，但我相信使用“Haversine 公式”之类的方法只会得到大约 0.5% 的平均误差，这在你的情况下可能是可接受的误差量。除非您谈论的是几英尺的距离，否则您总会有一些错误，即使在理论上地球也存在曲率......如果您需要更严格的 WGS-84 兼容方法，请查看“Vincenty 公式。 "

我了解starblue的来源，但好的软件工程通常需要权衡取舍，因此这完全取决于您所做工作所需的准确性。例如，从“曼哈顿距离公式”计算的结果与“距离公式”的结果相比，在某些情况下可能会更好，因为它的计算成本更低。想想“哪个点最近？” 不需要精确距离测量的场景。

关于，“Haversine 公式”很容易实现并且很好，因为它使用“球面三角”而不是基于“余弦定律”的基于二维三角的方法，因此您可以获得很好的准确性平衡过于复杂。

一位名叫Chris Veness的绅士有一个很棒的网站，它解释了您感兴趣的一些概念并演示了各种程序化实现；这也应该回答您的 x/y 转换问题。

于 2009-07-26T21:04:59.747 回答
13

在 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
于 2019-03-20T08:57:25.447 回答
9

转换GPS（WGS84)为笛卡尔坐标 的理论https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates

以下是我正在使用的：

  • GPS（WGS84）中的经度和笛卡尔坐标相同。
  • 纬度需要通过WGS 84椭球参数换算半长轴为6378137 m，和
  • 展平的倒数是 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高度。必须使用位势模型将MSL高度转换为h高于 的高度（Lemoine 等人，1998 年）。 这是通过对空间分辨率为 15 角分的大地水准面高度文件的网格进行插值来完成的。WGS 84 ellipsoidEGM96

或者，如果您有某个级别的专业人员 GPS具有海拔高度H（msl，高于平均海平面的高度）和，从内部表输出的所选数据的和UNDULATION之间的关系。你可以得到geoidellipsoid (m)h = H(msl) + undulation

通过笛卡尔坐标到 XYZ：

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

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

z = R *sin(lat)
于 2014-07-17T15:11:06.380 回答
6

proj.4软件提供了一个可以进行转换的命令行程序，例如

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

它还提供了一个C API。特别是，该函数pj_geodetic_to_geocentric将进行转换，而无需先设置投影对象。

于 2017-01-21T08:15:53.753 回答
4

如果您关心基于椭球而不是球体获取坐标，请查看Geographic_coordinate_conversion - 它提供了公式。GEodetic Datum具有转换所需的 WGS84 常数。

那里的公式还考虑了相对于参考椭球表面的高度（如果您从 GPS 设备获取高度数据，则很有用）。

于 2011-08-05T00:41:23.190 回答
3

为什么要实施已经实施和测试证明的东西？

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

具体来说，您的计算存在严重缺陷。地球不是一个完美的球体，地球半径的近似值可能无法精确测量。

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

于 2009-07-26T20:09:37.533 回答
1 
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);
于 2011-08-18T01:48:42.437 回答
1

我在 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])
于 2021-10-17T13:06:00.397 回答
0

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


}
于 2018-08-09T06:12:31.723 回答