python - 在两个经纬度点之间生成更多点

Question

在二维平面中，给定两个点(x1, y1)和，沿着两点之间的直线(x2, y2)生成等距点是直接的。N这也适用于 3D 平面。

但是，我正在尝试弄清楚您将如何处理地理坐标点。为了进一步说明我的观点，假设你有一个点 A(latA, lonA)代表它的纬度和经度，另一个点 B 代表它的纬度和经度(latB, lonB)。您将如何N在 A 和 B 之间生成点？是否有一个简单的库python可以实现这一目标？

标签： pythonpython-3.xlatitude-longitude

Answer 1

您可以直接使用 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_)
    ]).T

def 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 equator
A = 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. ]))