www说
您可以直接使用 numpy 执行此操作。这个想法是使用 3D 空间的标准插值公式,如A + d * (B - A). 像这样计算的点位于 A 和 B 之间的弦上,但可以投影回球体。为了在角度上有一个均匀的分布,我们需要从角度到弦上距离的映射,就像这里的图中这显示了均匀间隔角度的弦位置,并使用以下代码生成以检查正确性,因为所有角度和三角函数都很容易弄乱。def embed_latlon(lat, lon): """lat, lon -> 3d point""" lat_, lon_ = np.deg2rad(lat), np.deg2rad(lon) r = np.cos(lat_) return np.array([ r * np.cos(lon_), r * np.sin(lon_), np.sin(lat_) ]).Tdef project_latlon(x): """3d point -> (lat, lon)""" return ( np.rad2deg(np.arcsin(x[:, 2])), np.rad2deg(np.arctan2(x[:, 1], x[:, 0])) )def _great_circle_linspace_3d(x, y, n): """interpolate two points on the unit sphere""" # angle from scalar product alpha = np.arccos(x.dot(y)) # angle relative to mid point beta = alpha * np.linspace(-.5, .5, n) # distance of interpolated point to center of sphere r = np.cos(.5 * alpha) / np.cos(beta) # distance to mid line m = r * np.sin(beta) # interpolation on chord chord = 2. * np.sin(.5 * alpha) d = (m + np.sin(.5 * alpha)) / chord points = x[None, :] + (y - x)[None, :] * d[:, None] return points / np.sqrt(np.sum(points**2, axis=1, keepdims=True))def great_circle_linspace(lat1, lon1, lat2, lon2, n): """interpolate two points on the unit sphere""" x = embed_latlon(lat1, lon1) y = embed_latlon(lat2, lon2) return project_latlon(_great_circle_linspace_3d(x, y, n))# example on equatorA = 0, 0.B = 0., 30.great_circle_linspace(*A, *B, n=5)(array([0., 0., 0., 0., 0.]), array([ 0. , 7.5, 15. , 22.5, 30. ]))