如何在 WGS84 椭球上插值 ECEF 坐标答案

【问题标题】：How to interpolate ECEF coordinates on an WGS84 ellipsoid如何在 WGS84 椭球上插值 ECEF 坐标
【发布时间】：2014-02-10 20:01:38
【问题描述】：

是否有直接的方法（不涉及将坐标转换为纬度/经度）在 2 个 ECEF 坐标 (xyz) 之间进行插值，以便插值点位于 WGS84 椭球体上。原始的 2 个点是根据大地坐标计算得出的。

在球体上插值似乎很明显，但我似乎无法为椭球推导出解决方案。

提前谢谢你。

【问题讨论】：

沿什么插值？设置端点是不够的。你在插值什么函数？问题不完整。
我不确定我的问题中缺少什么；我试图在沿 WGS84 椭球的 2 个笛卡尔坐标之间进行插值。我找不到更清楚的方法。
@Philippe 使用将问题缩放为球形问题添加了答案，这极大地简化了事情（类似于将椭圆问题转换为圆形）。

标签： coordinates gis interpolation wgs84

【解决方案1】：

假设您有 2 个点 p0(x,y,z) 和 p1(x,y,z) 并希望在两者之间插入一些 p(t) 其中 t=<0.0,1.0>。

你可以：

将椭圆体重新缩放为球体

就像这样：
```
const double mz=6378137.00000/6356752.31414; // [m] equatoreal/polar radius of Earth
p0.z*=mz;
p1.z*=mz;
```
现在您获得了指向球形地球模型的笛卡尔坐标。
插值

简单的线性插值就可以了
```
p(t) = p0+(p1-p0)*t
```
但是你还需要对地球曲率进行归一化处理：
```
r0 = |p0|
r1 = |p1|

p(t) = p0+(p1-p0)*t
r(t) = r0+(r1-r0)*t

p(t)*=r/|p(t)|
```
其中|p0| 表示向量p0 的长度。
重新缩放回椭球体

除以相同的值
```
p(t).z/=mz
```

这既简单又便宜，但插值路径不会有线性时间尺度。

这里是 C++ 示例：

void XYZ_interpolate(double *pt,double *p0,double *p1,double t)
    {
    const double  mz=6378137.00000/6356752.31414;
    const double _mz=6356752.31414/6378137.00000;
    double p[3],r,r0,r1;
    // compute spherical radiuses of input points
    r0=sqrt((p0[0]*p0[0])+(p0[1]*p0[1])+(p0[2]*p0[2]*mz*mz));
    r1=sqrt((p1[0]*p1[0])+(p1[1]*p1[1])+(p1[2]*p1[2]*mz*mz));
    // linear interpolation
    r   = r0   +(r1   -r0   )*t;
    p[0]= p0[0]+(p1[0]-p0[0])*t;
    p[1]= p0[1]+(p1[1]-p0[1])*t;
    p[2]=(p0[2]+(p1[2]-p0[2])*t)*mz;
    // correct radius and rescale back
    r/=sqrt((p[0]*p[0])+(p[1]*p[1])+(p[2]*p[2]));
    pt[0]=p[0]*r;
    pt[1]=p[1]*r;
    pt[2]=p[2]*r*_mz;
    }

和预览：

黄色方块是使用的p0,p1笛卡尔坐标，白色曲线是插入路径t=<0.0,1.0> ...

【讨论】：