Unidata / MetPy

    above = broadcast_indicies(xp, minv2, ndim, axis)

    below = broadcast_indicies(xp, minv2 - 1, ndim, axis)

    # Calculate interpolation for each variable

    for i in range(len(args)):

        var = variables[i][below] + ((log_x_array - log_xp[below]) /

                                     (log_xp[above] - log_xp[below])) * (variables[i][above] -

                                                                         variables[i][below])

        # set to nan if extrap=set_nan

        if extrap == 'set_nan':

            var[minv == xp.shape[axis]] = np.nan

            var[minv == 0] = np.nan

        # Output to list

        ret.append(var)

    if len(ret) == 1:

        return ret[0]

    else:

        return ret

def broadcast_indicies(x, minv, ndim, axis):

    """Calculate index values to properly broadcast index array within data array.

    See usage in log_interp.

    """

    ret = []

    for dim in range(ndim):

        if dim == axis:

            ret.append(minv)

        else:

            broadcast_slice = [np.newaxis] * ndim

            broadcast_slice[dim] = slice(None)

    # Log x and xp

    log_x_array = np.log(x_array)

    log_xp = np.log(xp[sorter])

    # Calculate value above interpolated value

    minv = np.apply_along_axis(np.searchsorted, axis, xp[sorter], x[sort_x])

    minv2 = np.copy(minv)

    # If extrap is none and data is out of bounds, raise value error

    if ((np.max(minv) == xp.shape[axis]) or (np.min(minv) == 0)) and extrap == 'none':

        raise ValueError('Interpolation point out of data bounds encountered')

    # Warn if interpolated values are outside data bounds, will make these the values

    # at end of data range.

    if np.max(minv) == xp.shape[axis]:

        warnings.warn('Interpolation point out of data bounds encountered')

        minv2[minv == xp.shape[axis]] = xp.shape[axis] - 1

    if np.max(minv) == 0:

        warnings.warn('Interpolation point out of data bounds encountered')

        minv2[minv == 0] = 1

    # Create empty output list

    ret = []

    # Get indicies for broadcasting arrays

    above = broadcast_indicies(xp, minv2, ndim, axis)

    below = broadcast_indicies(xp, minv2 - 1, ndim, axis)

    # Calculate interpolation for each variable