python - 大圆到点的距离

Question

我有三点：A、B 和 C（WGS84 中的纬度、经度坐标）

我想计算 C 和 A、B 之间的最小路径之间的距离。（这必须适用于任何纬度、经度，所以我无法进行预测）

pyproj 或geographiclib.geodesic 可以为我做到这一点吗？

编辑：不重复，因为我想在 WGS84 中工作，而不是在球体上工作。正如我所说，这必须适用于任何经纬度。我知道这个问题没有简单的解决方案，这就是为什么我问是否有一种简单的方法可以使用 pyproj 或 geodesic 之类的库来做到这一点。

Answer 1

这是一种方法。它不是矢量化的，可能效率不高。我还没有测试所有可能的情况。这是一个图形示例：

#!/usr/bin/env python3
__author__  = 'hrmd'

from    geographiclib.geodesic  import  Geodesic
import  numpy                   as      np

def dist_from_great_circle_to_point(a, b, p):
    '''
    Calculates the shortest distance between a point P and a great
    circle passing through A and B. The method is:
        1. Find azimuth for line AP.
        2. Find azimuth for line AB.
        3. Consider point M which is the mirror image of point P
            in the great circle. Find azimuth AM from the first
            two azimuths.
        4. Find great circle between P and M. The point of interest
            is half-way along this line.

    Input
    a       [lon, lat] of first point defining the great circle.
    b       [lon, lat] of second point defining the great circle.
    p       [lon, lat] of query point.

    Output
    h       [lon, lat] on great circle closest to query point.
    d       Perpendicular distance between great circle and point.
    m       [lon, lat] of mirror image point.
    '''

    # Ellipsoid.
    ell     = Geodesic.WGS84

    # ab    Geodesic from A to B.
    # az_ab Azimuth of geodesic from A to B.
    ab      = ell.Inverse(a[1], a[0], b[1], b[0])
    az_ab   = ab['azi1']

    # ap    Geodesic from A to P.
    # az_ap Azimuth of geodesic from A to P.
    # d_p   Distance between A and P along geodesic.
    ap      = ell.Inverse(a[1], a[0], p[1], p[0])
    az_ap   = ap['azi1']
    d_ap    = ap['s12']

    # Point M is the image of P on the other side of line AB.
    # az_am Azimuth of geodesic from A to M.
    # am    Geodesic from A to M.
    # m     [lon, lat] of point M.
    az_am   = az_ab + (az_ab - az_ap)
    am      = ell.Direct(a[1], a[0], az_am, d_ap)
    m       = np.array([am['lon2'], am['lat2']])

    # Point H is the half-way between P and M along a great circle.
    # pm    Geodesic from P to M.
    # d     Half of the distance from P to M.
    # ph    Geodesic from P to H.
    # h     [lon, lat] of point H.
    pm      = ell.Inverse(p[1], p[0], m[1], m[0])
    d       = pm['s12']/2.0
    ph      = ell.Direct(p[1], p[0], pm['azi1'], d)
    h       = np.array([ph['lon2'], ph['lat2']])

    return h, d, m

def test():

    import  matplotlib.pyplot       as      plt
    from    mpl_toolkits.basemap    import  Basemap

    a   = np.array([ 80.0,  30.0])
    b   = np.array([100.0,  50.0])

    p   = np.array([ 70.0,  45.0])

    h, d, m = dist_from_great_circle_to_point(a, b, p)

    fig = plt.figure()
    ax  = plt.gca()
    bm  = Basemap(  projection  = 'ortho',
                    lat_0       = 35.0,
                    lon_0       = 85.0)
    bm.drawcoastlines(linewidth=0.25)

    points          = [a, b, p, m, h]
    point_names     = ['A', 'B', 'P', 'M', 'H']
    point_colors    = ['k', 'k', 'k', 'r', 'r']
    point_label_xy  = [(-10, -10), (10, 10), (-10, 5), (10, -5), (10, -3)]

    lines           = [[a, b], [p, m]]
    line_colors     = ['k', 'r']

    for point, color, name, label_xy in zip(
            points, point_colors, point_names, point_label_xy):

        bm.scatter( point[0], point[1],
                    color       = color,
                    marker      = '.',
                    latlon      = True)

        plt.annotate(   name,
                        xy          = bm(point[0], point[1]),
                        xycoords    = 'data',
                        xytext      = label_xy,
                        textcoords  = 'offset points',
                        color       = color)

    for line, color in zip(lines, line_colors):

        bm.drawgreatcircle(
                    line[0][0], line[0][1],
                    line[1][0], line[1][1],
                    color = color)

    plt.show()

if __name__ == '__main__':

    test()