Wsrtjtyk

In a 2D plane, given two points (x1, y1) and (x2, y2), it is straight forward to generate N equally spaced points along the straight line between the two points. This also applies for 3D plane.

However, I'm trying to work out how would you do so for geo-coordinated points. To illustrate my point further, say you have point A with (latA, lonA) that represents the lattitude and longitude of its, and another point B with (latB, lonB). How would you generate N points between A and B? Is there a straightforward library in python that could achieve this?

You can do this directly with numpy. The idea is to use the standard interpolation formula for 3D space, like A + d * (B - A). Points computed like this lie on the chord between A and B but can be projected back to the sphere.

In order to have a uniform distribution over angles, we need the mapping from angles to distances on the chord, like in the figure here

This shows chord locations for uniformly spaced angles and was generated with the code below to check correctness, since all the angles and trigonometric functions are easy to mess up.

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. ]))