interp.py

import numpy as np

def interp_coords(coords, tol):
    '''
    Interpolates extra coordinates when the euclidean distance between adjacent
    points are larger than given tolerance.
    The function tries to densify a list of coordinates based on a simple measure.
    If two adjacent points in the input coordinate list are further away from each
    other than the specified tolerance, the list will be densied between those two
    points. Roughly speaking, a point will be inserted every `tol` between the
    points.
    Inputs:
    -------
    coords:         2xN numpy array with geodetic coordinates
    tol:            Euclidean distance tolerance.
    Returns:
    --------
    A densified numpy array of coordinates. If tol < 0 the input array is returned
    unchanged.
    '''
    if tol < 0:
        return coords

    x, y = coords
    x = np.array(x)
    y = np.array(y)

    l = len(x)
    dsts = np.sqrt((x[0:-1]-x[1:l])**2 + (y[0:-1]-y[1:l])**2)
    I = dsts > tol

    offset = 0
    xy = np.ndarray(shape=(2, 0))
    for i, b in enumerate(I):
        xy = np.append(xy, (x[i], y[i]))
        if not b:
            continue

        n = round(dsts[i] / tol/2)
        x1 = np.linspace(x[i], x[i+1], n+2)
        y1 = np.linspace(y[i], y[i+1], n+2)
        xy = np.append(xy, list(zip(x1[1:], y1[1:])))

    # the last coordinate is not added in the loop above
    #xy = np.append(xy, (x[i+1], y[i+1])).reshape((-1,2)).tolist()
    xy = xy.reshape((-1, 2)).tolist()

    return xy

def interp_coords_new(coords, tol):
    '''
    Interpolates extra coordinates when the euclidean distance between adjacent
    points are larger than given tolerance.
    The function tries to densify a list of coordinates based on a simple measure.
    If two adjacent points in the input coordinate list are further away from each
    other than the specified tolerance, the list will be densied between those two
    points. Roughly speaking, a point will be inserted every `tol` between the
    points.
    Inputs:
    -------
    coords:         2xN numpy array with geodetic coordinates
    tol:            Euclidean distance tolerance.
    Returns:
    --------
    A densified numpy array of coordinates. If tol < 0 the input array is returned
    unchanged.
    '''
    if tol < 0:
        return coords

    x, y = coords
    x = np.array(x)
    y = np.array(y)

    dsts = np.hypot(np.diff(x), np.diff(y))
    # Points to the *beginning* of the segment.
    I = np.where(dsts > tol)[0]

    offset = 0
    xy = []
    for i in I:
        # Add points that are already below tolerance.
        xy.append((x[offset:i], y[offset:i]))

        # Interpolate between points above tolerance.
        n = np.ceil(dsts[i] / tol)
        x1 = np.linspace(x[i], x[i + 1], n + 1)
        y1 = np.linspace(y[i], y[i + 1], n + 1)
        xy.append((x1[:-1], y1[:-1]))

        offset = i + 1


    xy.append((x[offset:], y[offset:]))

    xy = np.concatenate(xy, axis=1)

    return xy.T


x = [0, 1, 2, 8, 9, 5, 10]
y = [0, 0, 0, 0, 0, 0, 0]

print(x)
print(np.array(interp_coords((x, y), 2))[:, 0])
print(interp_coords_new((x, y), 2)[:, 0])