Interpolate values and replace with NaNs within a long gap?

user_3015186 0
Viewed 451

I am trying to interpolate data with gaps. Sometimes the gap can be very large, and I do not want the interpolation to "succeed" within the gap; the result should be NaNs inside a large gap. For example, consider this example data set:

orig_x = [26219, 26225, 26232, 28521, 28538]
orig_y = [39, 40, 41, 72, 71]

which has clear gap between x-values 26232 and 28521. Now, I would like to have the orig_y interpolated to x-values like this:

import numpy as np 

x_target = np.array(range(min(orig_x) // 10 * 10 + 10, max(orig_x) // 10 * 10 + 10, 10))
# array([26220, 26230, 26240, 26250, 26260, 26270, 26280, 26290,
#        ...
#       28460, 28470, 28480, 28490, 28500, 28510, 28520, 28530])

and the output y_target should be np.nan everywhere else than at 26220, 26230 and 28520. Let's say that the condition for this would be that if there is a gap larger than 40 in the data, the interpolation should result to np.nan inside this data gap.

Goal shown as a picture

Instead of this

not this

Get something like this

but this

i.e. the "gap" in the data should result to np.nan instead of garbage data.

Question

What would be the best way (fastest interpolation) to achieve this kind of interpolation? The interpolation can be linear or more sophisticated (e.g. cubic spline). One possibility I have in mind would be to use the scipy.interpolate.interp1d as starting point like this

from scipy.interpolate import interp1d
f = interp1d(orig_x, orig_y, bounds_error=False)
y_target = f(x_target)

and then search for gaps in the data and replace the interpolated data with np.nan inside the gaps. Since I will be using this on fairly large dataset (~10M rows, few columns, handled in parts), performance is a key.

interpolation python
1 Answers

After some trial and error, a think I got a "fast enough" implementation using basic linear interpolation and numba for speedups. Forgive for writing everything in the same loop and same function, but it seems that is the numba way of making your code fast. (numba loves loops, and does not seem to accept nested functions)

Test data used

I added some mode data to x_target to test the algorithm performance.

orig_x = np.array([26219, 26225, 26232, 28521, 28538])
orig_y = np.array([39, 40, 41, 72, 71])

x_target = np.array(
    np.arange(min(orig_x) // 10 * 10,
              max(orig_x) // 10 * 10 + 10, 0.1))

Test code

from matplotlib import pyplot as plt

y_target = interpolate_with_max_gap(orig_x, orig_y, x_target, max_gap=40)

plt.scatter(x_target, y_target, label='interpolated', s=10)
plt.scatter(orig_x, orig_y, label='orig', s=10)
plt.legend()
plt.show()

Test results

The data is interpolated in regions with gap less than max_gap (40): test results

closeup:

closeup

Speed:
I first tried a pure python + numpy implementation, which took 49.6 ms with the same test data (using timeit). This implementation with numba takes 480µs (100x speedup!). When using target_x_is_sorted=True, the speed is 80.1µs!

The orig_x_sorted=True did not give speedup, probably since the orig_x is so short that sorting it does not make any difference in timing in this example.

Implementation

import numba
import numpy as np

@numba.njit()
def interpolate_with_max_gap(orig_x,
                             orig_y,
                             target_x,
                             max_gap=np.inf,
                             orig_x_is_sorted=False,
                             target_x_is_sorted=False):
    """
    Interpolate data linearly with maximum gap. If there is
    larger gap in data than `max_gap`, the gap will be filled
    with np.nan.

    The input values should not contain NaNs.

    Parameters
    ---------
    orig_x: np.array
        The input x-data
    orig_y: np.array
        The input y-data
    target_x: np.array
        The output x-data; the data points in x-axis that
        you want the interpolation results from.
    max_gap: float
        The maximum allowable gap in `orig_x` inside which
        interpolation is still performed. Gaps larger than
        this will be filled with np.nan in the output `target_y`.
    orig_x_is_sorted: boolean, default: False
        If True, the input data `orig_x` is assumed to be monotonically
        increasing. Some performance gain if you supply sorted input data.
    target_x_is_sorted: boolean, default: False
        If True, the input data `target_x` is assumed to be 
        monotonically increasing. Some performance gain if you supply
        sorted input data.

    Returns
    ------
    target_y: np.array
        The interpolation results.
    """
    if not orig_x_is_sorted:
        # Sort to be monotonous wrt. input x-variable.
        idx = orig_x.argsort()
        orig_x = orig_x[idx]
        orig_y = orig_y[idx]

    if not target_x_is_sorted:
        target_idx = target_x.argsort()
        # Needed for sorting back the data.
        target_idx_for_reverse = target_idx.argsort()
        target_x = target_x[target_idx]

    target_y = np.empty(target_x.size)
    idx_orig = 0
    orig_gone_through = False

    for idx_target, x_new in enumerate(target_x):

        # Grow idx_orig if needed.
        while not orig_gone_through:

            if idx_orig + 1 >= len(orig_x):
                # Already consumed the orig_x; no more data
                # so we would need to extrapolate
                orig_gone_through = True
            elif x_new > orig_x[idx_orig + 1]:
                idx_orig += 1
            else:
                # x_new <= x2
                break

        if orig_gone_through:
            target_y[idx_target] = np.nan
            continue

        x1 = orig_x[idx_orig]
        y1 = orig_y[idx_orig]
        x2 = orig_x[idx_orig + 1]
        y2 = orig_y[idx_orig + 1]

        if x_new < x1:
            # would need to extrapolate to left
            target_y[idx_target] = np.nan
            continue

        delta_x = x2 - x1

        if delta_x > max_gap:
            target_y[idx_target] = np.nan
            continue

        delta_y = y2 - y1

        if delta_x == 0:
            target_y[idx_target] = np.nan
            continue

        k = delta_y / delta_x

        delta_x_new = x_new - x1
        delta_y_new = k * delta_x_new
        y_new = y1 + delta_y_new

        target_y[idx_target] = y_new

    if not target_x_is_sorted:
        return target_y[target_idx_for_reverse]
    return target_y
user_30151860
Related