有了解 GIS 的么，请教个经纬度的问题？

[James369]

在地球上给定 2 个点 a 和 b 的经纬度，沿着球面大圆连接 ab 两个点做一条最短线段。
现在求解：在 ab 上等距离的划分 n 个点，求出这 n 个点的经纬度坐标？

还有 2 个细节问题：
1. 地球不是一个真正的球体，它是南北两极稍扁的椭球体（半径差了几十公里）。
2. ab 连线理论上可能通过极点（极点在某些空间 /投影变换的时候可能消失），那么软件层面应该如何处理这种极端情况。

[anytk]

参考： https://geographiclib.sourceforge.io/html/python/geodesics.html#introduction
先 Geodesic.Inverse 然后距离平分再 Geodesic.Direct

[liyang5945]

之前用过 turf.js ，需要用里面三个方法，先算出 ab 两点的距离，比如 100Km，然后计算 ab 之间的方位角，比如 45 度(即 b 在 a 的东北方向，正北为 0 度)，再根据点、距离和角度计算目标点(以 a 为基准点，45 度上任意距离点的坐标)
文档地址： https://turfjs.fenxianglu.cn/category/measurement/destination.html

[VensonEEE]

先切经线 再切纬线

[ck65]

看看 linear referencing 的知识，你的问题不复杂。比如 PostGIS 的 line interpolate points 函数或许可以直接解决你的问题 https://postgis.net/docs/ST_LineInterpolatePoints.html

[urmyfaith]

先转换空间，然后等分，再将空间转换回来。

[James369]

@anytk 有一点点思路，就是把经纬度坐标转化为 球坐标（类似弧度）处理。

[sadfQED2]

先把经纬度坐标系转换成墨卡托投影坐标系，然后再进行几何计算，计算完以后再转换成经纬度坐标

[3dwelcome]

不用那么复杂，四元数就是专门干在球面上，两点之间插值算法的。线性插值就是了。

[no1xsyzy]

用坐标系变换做会更简单一点。

原点不变，不作切变，只做旋转，三个自由度
a 点变换到 $(1,0,0)$
b 点变换到 $(\cos\theta, \sin\theta, 0)$
4 个变量四个方程（实际上是 6 个方程但有两个同解）
更直接的求解是按 x=y=0 轴把 a 转到 y=0 平面，再按 x=z=0 轴把 a 转到 z=0 平面，再用 y=z=0 轴把 b 转到 z=0 平面。
从而得知等分后的点为 $(\cos\frac{k\theta}{n+1}}, \sin\frac{k\theta}{n+1}}, 0)$
之后再旋转回去。

[3dwelcome]

@sadfQED2 墨卡托平面投影里，投影后的经纬度网格是曲线。
如果平面里你画一条直线，那反投影到球面上，也是弯的，肯定不是球面上的最短距离。
只有四元数才是最短距离。

[no1xsyzy]

@sadfQED2 墨卡托投影不能平分距离。
@3dwelcome 确实三维空间内四元数比矩阵有不少优势。
@anytk 这个是椭球模型吧，比楼主的问题更一般化。

[James369]

@3dwelcome 怎么个插值法呢，

[no1xsyzy]

@3dwelcome 墨卡托没有破坏直线，方法是投影到赤道圈沿极轴无限延伸形成的圆柱，然后沿任意经线切开或者循环表示。
但是破坏了距离，所以不能用于等分插值。简单的示例是极点到赤道上任意一点的二等分插值，应在 45 度纬度，但墨卡托下极点在无穷远处，平面插值二分之一无穷远也是无穷远。

[no1xsyzy]

@James369 en.wikipedia.org/wiki/Slerp

[KouShuiYu]

先将经纬度换算成 xyz 坐标，然后当球体等分处理就行了

[ChristianChen]

两点分别算经纬度弧长，平分后再反求经差，纬差；投影那是另外一个问题
公式可以参阅 http://www.whigg.cas.cn/resource/zysk/dlxx/201011/P020101127664602088413.pdf
P15:2-13，P16:2-17
自己验证吧，我已经完全没有记忆了

[cht]

写了个简单 demo，如果有 Xcode 的话可以试一下：aHR0cHM6Ly9naXRodWIuY29tL3lvMTk5NS9EYWlseV9Td2lmdF9UYXNrcy90cmVlL21hc3Rlci9HZW9kZXRpY0RpdmlkZUFMaW5l

仅供参考