Быстрая 1D линейная np.NaN интерполяция над большой 3D-матрицей

Question

У меня есть 3D-массив (z, y, x) с shape=(92, 4800, 4800), где каждое значение вдоль axis 0 представляет собой другой момент времени. Приобретение значений во временном домене не удалось в нескольких случаях, вызвав некоторые значения np.NaN. В других случаях значения не были получены, и все значения вдоль z составляют np.NaN.

Каков наиболее эффективный способ использования линейной интерполяции для заполнения np.NaN по адресу axis 0 без учета экземпляров, где все значения равны np.NaN?

Вот рабочий пример того, что я делаю, который использует pandas wrapper для scipy.interpolate.interp1d. Это занимает около 2 секунд на каждый фрагмент исходного набора данных, что означает, что весь массив обрабатывается в 2.6 часа. Примерный набор данных с уменьшенным размером занимает около 9,5 секунд.

import numpy as np 
import pandas as pd 

# create example data, original is (92, 4800, 4800) 
test_arr = np.random.randint(low=-10000, high=10000, size=(92, 480, 480)) 
test_arr[1:90:7, :, :] = -32768 # NaN fill value in original data 
test_arr[:, 1:90:6, 1:90:8] = -32768 

def interpolate_nan(arr, method="linear", limit=3): 
    """return array interpolated along time-axis to fill missing values""" 
    result = np.zeros_like(arr, dtype=np.int16) 

    for i in range(arr.shape[1]): 
     # slice along y axis, interpolate with pandas wrapper to interp1d 
     line_stack = pd.DataFrame(data=arr[:,i,:], dtype=np.float32) 
     line_stack.replace(to_replace=-37268, value=np.NaN, inplace=True) 
     line_stack.interpolate(method=method, axis=0, inplace=True, limit=limit) 
     line_stack.replace(to_replace=np.NaN, value=-37268, inplace=True) 
     result[:, i, :] = line_stack.values.astype(np.int16) 
    return result 

Производительность на моей машине с набором данных, например:

%timeit interpolate_nan(test_arr) 
1 loops, best of 3: 9.51 s per loop 

Edit:

я должен уточнить, что код производит мой ожидаемый результат. Вопрос в том, как я могу оптимизировать этот процесс?

Answer 1

Я недавно решил эту проблему для своего конкретного случая использования с помощью numba, а также сделал a little writeup on it.

from numba import jit 

@jit(nopython=True) 
def interpolate_numba(arr, no_data=-32768): 
    """return array interpolated along time-axis to fill missing values""" 
    result = np.zeros_like(arr, dtype=np.int16) 

    for x in range(arr.shape[2]): 
     # slice along x axis 
     for y in range(arr.shape[1]): 
      # slice along y axis 
      for z in range(arr.shape[0]): 
       value = arr[z,y,x] 
       if z == 0: # don't interpolate first value 
        new_value = value 
       elif z == len(arr[:,0,0])-1: # don't interpolate last value 
        new_value = value 

       elif value == no_data: # interpolate 

        left = arr[z-1,y,x] 
        right = arr[z+1,y,x] 
        # look for valid neighbours 
        if left != no_data and right != no_data: # left and right are valid 
         new_value = (left + right)/2 

        elif left == no_data and z == 1: # boundary condition left 
         new_value = value 
        elif right == no_data and z == len(arr[:,0,0])-2: # boundary condition right 
         new_value = value 

        elif left == no_data and right != no_data: # take second neighbour to the left 
         more_left = arr[z-2,y,x] 
         if more_left == no_data: 
          new_value = value 
         else: 
          new_value = (more_left + right)/2 

        elif left != no_data and right == no_data: # take second neighbour to the right 
         more_right = arr[z+2,y,x] 
         if more_right == no_data: 
          new_value = value 
         else: 
          new_value = (more_right + left)/2 

        elif left == no_data and right == no_data: # take second neighbour on both sides 
         more_left = arr[z-2,y,x] 
         more_right = arr[z+2,y,x] 
         if more_left != no_data and more_right != no_data: 
          new_value = (more_left + more_right)/2 
         else: 
          new_value = value 
        else: 
         new_value = value 
       else: 
        new_value = value 
       result[z,y,x] = int(new_value) 
    return result 

Это примерно 20 раз быстрее, чем мой исходного кода.

Answer 2

Это зависит от многого. вам придется вытащить лист бумаги и рассчитать погрешность вашей общей статистики, если вы не интерполировать и просто заполнить нулями NaN.

Кроме этого, я думаю, что ваша интерполяция находится сверху. Просто найдите каждый NaN и линейно интерполируйте соседние четыре значения (что суммирует значения в (y + - 1, x + - 1)) - это серьезно ограничит вашу ошибку (вычислите себя!), и вы не интерполируете с каким-либо сложным методом в вашем случае (вы не определили method).

Вы можете попытаться просто предварительно вычислить одну «усредненную» матрицу 4800x4800 за значение z - это не должно занять много времени - путем применения кросс-формы ядра по всей матрице (все это очень обработка изображений - как здесь). В случае NaN некоторые из усредненных значений будут NaN (каждый усредненный пиксель, где NaN был в соседстве), но вам все равно - если нет двух соседних NaN, ячеек NaN, которые вы хотите заменить в исходная матрица все вещественнозначна.

Затем вы просто заменяете все NaN на значение в усредненной матрице.

Сравните скорость этого со скоростью «ручного» расчета среднего значения по окрестностям для каждого найденного вами NaN.