从经纬度转换为笛卡尔坐标系

Question

我有一些以纬度和经度（WGS-84）表示的以地球为中心的坐标点。

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

- daphshez

---

你有成功将 WGS-84 经纬度转换为笛卡尔坐标系吗？我还有海拔高度信息。我已经尝试了这里被接受的答案，但是它给出的结果并不正确。我和这个网站 https://www.apsalin.com/convert-geodetic-to-cartesian.aspx 比较了一下我的结果。 - Yasmin

Answer 1

这是我找到的答案：

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

- x轴经过经纬度(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#代码（请注意，它与公式非常不同），还有一些解释和漂亮的图表，说明为什么这是正确的。

- daphshez

---

错误了，你假设地球是一个球体，而WGS-84假设它是一个椭球体。 - starblue

@starblue: 我不确定您有没有资格将给定答案标记为“正确”或“错误”。使用可用的参考WGS-84经纬度，进行球面逼近（以获取ECEF样式的x、y、z坐标）要么对原贴提问者的需求“足够”，要么“不足够”。对于距离和方位角估计，我觉得这种简单转换足够了。如果他正在发射卫星，可能就不行了。毕竟，WGS-84本身是“错误的”...因为它不是地球表面的完美模型；所有椭球模型都是近似值。很遗憾原帖提问者没有告诉我们他打算做什么。 - Dan H

问题要求使用WGS-84，如果您回答其他内容，您至少应讨论差异/误差，而这个回答没有做到。 - starblue

@daphna-shezaf 无法进行反向转换...我已经将弧度转换为角度，但结果不同... - user2402179

不使用 z，只有起始点的 x 和 y，是否可能进行从笛卡尔坐标到常规坐标的转换？ - osk

Answer 2

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

是的，WGS-84假设地球是一个椭球体，但我认为使用“Haversine公式”这样的方法平均误差只有大约0.5%，在你的情况下可能是可接受的误差范围。除非你要测量几英尺的距离，否则总会存在一些误差，即使在那种情况下理论上也有地球曲率... 如果你需要更严格符合WGS-84标准的方法，请查看“Vincenty公式”。

我理解starblue的观点，但优秀的软件工程通常涉及权衡，因此它完全取决于你所需的精度。例如，“曼哈顿距离公式”和“距离公式”计算出来的结果在某些情况下可以更好，因为计算成本更低。想象一下需要“哪个点最近？”这种情况，在这种情况下你不需要精确的距离测量。

关于“Haversine公式”，它易于实现，并且很好，因为它使用了“球面三角法”，而不是基于二维三角法的“余弦定理”方法，因此你可以在精度和复杂性之间获得一个很好的平衡。

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

- bn.

---

0.5%的误差 - 0.5%是什么？在这个问题的背景下，它可能是地球的半径，因此0.5%可能是30公里 :) - MarkJ

查看了你的链接。0.5% 的引用是指两点之间大圆距离误差，因此与这个问题没有直接关系。我认为，在将纬度和经度转换为以地球中心为原点的笛卡尔坐标时，假设地球是球形的误差可能会很大。不清楚提问者想要用笛卡尔坐标做什么。要么出于某种奇怪的原因，使用笛卡尔坐标更方便，要么可能是数据导出的要求？如果是后者，准确性就很重要。 - MarkJ

Answer 3

在Python 3.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

- Mayank Kumar

---

如果lat，lon的值是以十进制表示（即以米为单位），我们需要将R值从千米转换为米来计算x和y吗？ - user2293224

我使用这种方法在CesiumJS中绘制点到地球上，它使用笛卡尔坐标而不是经纬度/高度。我还需要将公里转换为米...很简单！只需在使用x、y和z变量的每行末尾添加* 1000即可。 - Union In Design

Answer 4

将GPS（WGS84）转换为笛卡尔坐标的理论 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椭球体的高度。

通常，GPS会给出相对于MSL高度的H。 必须使用地势势能模型EGM96（Lemoine et al, 1998）将MSL高度转换为相对于WGS 84椭球体的高度h。
这是通过插值一个空间分辨率为15角分的大地水准面高度文件的网格来完成的。

或者，如果您有一些专业级别的GPS，具有高度H（msl，相对于平均海平面的高度）和UNDULATION，则可以通过内部表格中所选数据输出的大地水准面和椭球体（m）之间的关系来获得h = H(msl) + undulation 使用笛卡尔坐标转换为XYZ坐标：

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

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

z = R *sin(lat)

- Howie

---

R的值是多少？ - eych

我猜想这是球体的半径，对于地球而言它是6371公里。 - Matthias

Answer 5

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

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

它还提供了一个C API。上述的C等效代码可能是这样的

#include <math.h>
#include <proj.h>
#include <stdio.h>

int main(int argc, char* argv[])
{

    PJ* P = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
        "+proj=latlong +datum=WGS84",
        "+proj=geocent +datum=WGS84", NULL);

    if (P == NULL) {
        return 1;
    }

    PJ_COORD c;
    c.lpzt.lam = -110.0;  /* Longitude in degrees */
    c.lpzt.phi = 40.0;    /* Latitude in degrees */
    c.lpzt.z = 0.0;       /* Height in meters */
    c.lpzt.t = HUGE_VAL;  /* Time (not used here) */

    PJ_COORD c_out = proj_trans(P, PJ_FWD, c);

    printf("X = %.17g m\n", c_out.xyzt.x);
    printf("Y = %.17g m\n", c_out.xyzt.y);
    printf("Z = %.17g m\n", c_out.xyzt.z);

    return 0;
}

导致输出结果

X = -1673404.5546274509 m
Y = -4597641.2274514409 m
Z = 4077985.5722003761 m

还有一个可用的Python包：

from pyproj import Transformer
lon, lat, hgt = -110., 40., 0.
llh_to_xyz = Transformer.from_crs("+proj=latlong +datum=WGS84", "+proj=geocent +datum=WGS84")
print(llh_to_xyz.transform(lon, lat, hgt))

- Brian Hawkins

---

看起来pj_geodetic_to_geocentric()不再在API中。 - Chris Gnam

似乎pj_geodetic_to_geocentric在2019年随着PROJ 6版本的发布被弃用了。我已经更新了答案，提供了一个适用于PROJ 9.2.0版本的C示例。 - Brian Hawkins

似乎pj_geodetic_to_geocentric在2019年随着PROJ 6版本的发布被弃用了。我已经更新了答案，提供了一个适用于PROJ 9.2.0版本的C示例。 - undefined

Answer 6

如果您关心的是基于椭球体而不是球体获取坐标，请参考地理坐标转换 - 它提供公式。大地基准具有进行转换所需的WGS84常数。

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

- Stjepan Rajko

Answer 7

我在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])

- Tom Lakeside

---

你能具体说明一下 h 参数是什么吗？它是指点相对于海平面的高度还是到地球中心的距离？ - Krystof

它是海拔高度。 - Tom Lakeside

Answer 8

在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;


}

- Ashutosh Chapagain

---

参数 alt 是什么？ - baliman

高度，如果你不知道GPS是如何工作的，你在这里干嘛呢 ;) - MushyPeas

问题中根本没有提到GPS或高度。 - teekarna

Answer 9

为什么要实现已经实现和经过测试的东西呢？

比如说，C#有NetTopologySuite，它是JTS Topology Suite的.NET端口。

具体来说，你的计算中存在严重缺陷。地球并不是一个完美的球体，对地球半径的近似可能无法满足精确测量的要求。

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

- Yuval Adam

---

+1. 使用可靠的库比自行编写函数更准确且更容易。 - MarkJ

NetTopologySuite 如何将经纬度转换成笛卡尔坐标系？ - vinayan

NTS不包括坐标转换功能，也许你需要使用Proj.NET http://projnet.codeplex.com - D_Guidi

荒谬，答案甚至没有提供转换能力。 - Motomotes

Answer 10

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

- yang jun

---

你能详细说明一下吗？我创建了一个简单的应用程序，试图使用你的方法转换单个坐标。但它总是失败，因为源（2）的维度和目标（3）的维度不同，导致异常java.lang.IllegalArgumentException: dimension must be <= 3。 - oschrenk

嗯...我稍微看了一下JTS。直到及包括新的LineString()那一行看起来像是JTS。但我没有在JTS中看到CRS和Transform相关的内容。所以：它们在那里而我却没有发现吗？它们曾经存在且在1.12中被移除了吗？还是说这是另一个库的功能？ - Dan H