Линейная интерполяция Python 4D на прямоугольной сетке

Мне нужно интерполировать данные температуры линейно в 4 измерениях (широта, долгота, высота и время).
Количество точек довольно высокое (360x720x50x8), и мне нужен быстрый метод вычисления температуры в любой точке пространства и времени в пределах границ данных.

Я попытался использовать scipy.interpolate.LinearNDInterpolator, но использование Qhull для триангуляции неэффективно на прямоугольной сетке и занимает несколько часов.

Прочитав этот билет SciPy, решение, похоже, реализовало новый nd-интерполятор с использованием стандартного interp1d для вычисления большего числа точек данных, а затем использовать подход "ближайшего соседа" с новым набором данных.

Это, однако, занимает много времени (минуты).

Есть ли быстрый способ интерполяции данных по прямоугольной сетке в 4 измерениях без выполнения минут?

Я думал использовать interp1d 4 раза, не вычисляя более высокую плотность точек, но оставляя его для вызова пользователем с координатами, но я не могу понять, как это сделать.

В противном случае я бы написал здесь свой собственный 4D-интерполятор, специфичный для моих потребностей?

Вот код, который я использовал для тестирования:

Использование scipy.interpolate.LinearNDInterpolator:

import numpy as np
from scipy.interpolate import LinearNDInterpolator

lats = np.arange(-90,90.5,0.5)
lons = np.arange(-180,180,0.5)
alts = np.arange(1,1000,21.717)
time = np.arange(8)
data = np.random.rand(len(lats)*len(lons)*len(alts)*len(time)).reshape((len(lats),len(lons),len(alts),len(time)))

coords = np.zeros((len(lats),len(lons),len(alts),len(time),4))
coords[...,0] = lats.reshape((len(lats),1,1,1))
coords[...,1] = lons.reshape((1,len(lons),1,1))
coords[...,2] = alts.reshape((1,1,len(alts),1))
coords[...,3] = time.reshape((1,1,1,len(time)))
coords = coords.reshape((data.size,4))

interpolatedData = LinearNDInterpolator(coords,data)

Используя scipy.interpolate.interp1d:

import numpy as np
from scipy.interpolate import LinearNDInterpolator

lats = np.arange(-90,90.5,0.5)
lons = np.arange(-180,180,0.5)
alts = np.arange(1,1000,21.717)
time = np.arange(8)
data = np.random.rand(len(lats)*len(lons)*len(alts)*len(time)).reshape((len(lats),len(lons),len(alts),len(time)))

interpolatedData = np.array([None, None, None, None])
interpolatedData[0] = interp1d(lats,data,axis=0)
interpolatedData[1] = interp1d(lons,data,axis=1)
interpolatedData[2] = interp1d(alts,data,axis=2)
interpolatedData[3] = interp1d(time,data,axis=3)

Большое спасибо за вашу помощь!

Ответ 1

В том же билете, который вы связали, есть пример реализации того, что они называют интерполяцией тензорного продукта, показывая правильный способ вложения рекурсивных вызовов в interp1d. Это эквивалентно квадрилинейной интерполяции, если вы выбираете параметр kind='linear' по умолчанию для interp1d.

Хотя это может быть достаточно хорошим, это не линейная интерполяция, и в интерполяционной функции будут более высокие порядковые члены, так как это изображение из записи в википедии на билинейной интерполяции показывает:

Это может быть очень хорошо для того, что вам нужно, но есть приложения, в которых предпочтительна триангулированная, действительно кусочно-линейная интерполяция. Если вам это действительно нужно, есть простой способ работать вокруг медленности qhull.

Как только LinearNDInterpolator был установлен, есть два шага, чтобы приблизиться к интерполированному значению для данной точки:

выяснить внутри, в каком треугольнике (4D гипертетраэдр в вашем случае) точка есть, и
интерполировать с помощью барицентрических координат точки относительно вершин как веса.

Вероятно, вы не хотите связываться с барицентрическими координатами, поэтому лучше оставить это на LinearNDInterpolator. Но вы знаете кое-что о триангуляции. В основном, поскольку у вас есть регулярная сетка, внутри каждого гиперкуба триангуляция будет одинаковой. Таким образом, чтобы интерполировать одно значение, вы можете сначала определить, в каком субкубе ваша точка, постройте LinearNDInterpolator с 16 вершинами этого куба и используйте его для интерполяции вашего значения:

from itertools import product

def interpolator(coords, data, point) :
    dims = len(point)
    indices = []
    sub_coords = []
    for j in xrange(dims) :
        idx = np.digitize([point[j]], coords[j])[0]
        indices += [[idx - 1, idx]]
        sub_coords += [coords[j][indices[-1]]]
    indices = np.array([j for j in product(*indices)])
    sub_coords = np.array([j for j in product(*sub_coords)])
    sub_data = data[list(np.swapaxes(indices, 0, 1))]
    li = LinearNDInterpolator(sub_coords, sub_data)
    return li([point])[0]

>>> point = np.array([12.3,-4.2, 500.5, 2.5])
>>> interpolator((lats, lons, alts, time), data, point)
0.386082399091

Это не может работать с векторизованными данными, поскольку для этого требуется хранить LinearNDInterpolator для каждого возможного субкуба, и хотя он, вероятно, будет быстрее, чем триангулировать все это, он все равно будет очень медленным.

Ответ 2

scipy.ndimage.map_coordinates является хорошим быстрым интерполятором для однородных сеток (все ящики одинакового размера). См. multivariate-spline-interpolation-in-python-scipy на SO для четкого описания.

Для неравномерных прямоугольных сеток простая обертка Intergrid отображает/масштабирует неравномерные равномерные сетки, то map_координаты. На 4d тестовом примере, таком как ваш, для каждого запроса требуется около 1 & mu; sec:

Intergrid: 1000000 points in a (361, 720, 47, 8) grid took 652 msec

Ответ 3

Для очень похожих вещей я использую Scientific.Functions.Interpolation.InterpolatingFunction.

    import numpy as np
    from Scientific.Functions.Interpolation import InterpolatingFunction

    lats = np.arange(-90,90.5,0.5)
    lons = np.arange(-180,180,0.5)
    alts = np.arange(1,1000,21.717)
    time = np.arange(8)
    data = np.random.rand(len(lats)*len(lons)*len(alts)*len(time)).reshape((len(lats),len(lons),len(alts),len(time)))

    axes = (lats, lons, alts, time)
    f = InterpolatingFunction(axes, data)

Теперь вы можете оставить его пользователю, чтобы вызвать InterpolatingFunction с координатами:

>>> f(0,0,10,3)
0.7085675631375401

InterpolatingFunction имеет приятные дополнительные функции, такие как интеграция и нарезка.

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