Интерполяция 3d-массива в Python расширена

Я знаю, как вы можете избавиться от этих циклов, но вам это не понравится.

Проблема заключается в том, что это использование interp1d дает по существу матрично- interp1d функцию, интерполированную в 1d-домене, то есть функцию F(x) где для каждого скаляра x вас есть 2-мерная форма F Интерполяция, которую вы пытаетесь сделать, заключается в следующем: создание индивидуального интерполятора для каждой из вас (lat,lon) . Это больше по линиям F_(lat,lon)(x) .

Причина этого в том, что для вашего случая использования вы вычисляете матрицу F(x) для каждой из ваших точек запроса, но затем продолжаете отбрасывать все элементы матрицы, за исключением одного (элемент [lat,lon] для точки запроса, соответствующей этой паре). Таким образом, вы делаете кучу ненужных вычислений, вычисляя все эти нерелевантные значения функций. Проблема в том, что я не уверен, что есть более эффективный способ.

Ваш случай использования может быть исправлен с соответствующей памятью за спиной. Тот факт, что ваши циклы работают часами, говорит о том, что это будет невозможно для вашего варианта использования, но в любом случае я покажу это решение. Идея состоит в том, чтобы превратить ваш 3D-массив в 2-й, сделать интерполяцию с этой формой, а затем взять диагональные элементы вдоль эффективного 2d-подпространства вашего интерполированного результата. Вы все равно будете вычислять каждый нерелевантный матричный элемент для каждой точки запроса, но вместо того, чтобы перебирать ваши массивы, вы сможете извлечь соответствующие элементы матрицы с помощью одного шага индексации. Стоимость этого заключается в создании гораздо большего вспомогательного массива, который, скорее всего, не поместится в вашу свободную оперативную память.

Во всяком случае, вот трюк в действии, сравнивающий ваш текущий подход с одним:

 import numpy as np from scipy.interpolate import interp1d arr = np.random.randint(0, 9, size=(3, 4, 5)) x = np.linspace(0, 10, 3) x_new = np.random.rand(6,4,5)*10 ## x is shape 3 ## arr is shape 3x4x5 ## x_new is shape 6x4x5 # original, loopy approach interp_func = interp1d(x, arr, axis=0) res = np.empty((6, 4, 5)) for lat in range(res.shape[1]): for lon in range(res.shape[2]): temp = interp_func(x_new[:,lat,lon]) # shape (6,4,5) each iteration res[:,lat,lon] = temp[:,lat,lon] # new, vectorized approach arr2 = arr.reshape(arr.shape[0],-1) # shape (3,20) interp_func2 = interp1d(x,arr2,axis=0) x_new2 = x_new.reshape(x_new.shape[0],-1) # shape (6,20) temp = interp_func2(x_new2) # shape (6,20,20): 20 larger than original! s = x_new2.shape[1] # 20, used for fancy indexing ranges res2 = temp[:,range(s),range(s)].reshape(res.shape) # shape (6,20) -> (6,4,5) 

Результирующие массивы res и res2 абсолютно одинаковы, поэтому подход, вероятно, работает. Но, как я уже сказал, для массива запросов (nx,nlat,nlon) нам нужен вспомогательный массив формы (nx,nlat*nlon,nlat*nlon) , который, как правило, будет иметь огромную потребность в памяти.

Единственная строгая альтернатива, которую я могу придумать, – это просто выполнить ваши одномерные интерполяции один за другим: определение интерполяторов nlat*nlon в двойном цикле. Это будет иметь гораздо большие накладные расходы на создание интерполяторов, но, с другой стороны, вы не будете делать кучу лишних рабочих вычислений интерполированных значений массива, которые вы затем отбросите.

Наконец, в зависимости от вашего использования ase вы должны рассмотреть возможность использования многомерной интерполяции (я думаю, interpolate.interpnd или interpolate.griddata ). Предполагая, что ваша функция является гладкой как функция широты и долготы, возможно, имеет смысл интерполировать ваш полный набор данных в более высоком измерении. Таким образом, вам нужно создать свой интерполятор один раз и запросить именно то, что вам нужно, без лишнего пуха на вашем пути.

Если вы в конечном итоге придерживаетесь своей текущей реализации, вы, вероятно, можете значительно повысить производительность, перемещая ось интерполяции до последней позиции. Таким образом, всякая векторная операция действует на смежных блоках памяти (предполагая порядок памяти по умолчанию C), и это хорошо согласуется с философией «сборка 1-го массива». Поэтому вы должны что-то сделать в соответствии с

 arr = arr.transpose(1,2,0) # shape (4,5,3) interp_func = interp1d(x, arr, axis=-1) ... for lat ...: for lon ...: res[lat,lon,:] = temp[lat,lon,:] # shape (4,5,6) 

Если вам нужно восстановить исходный порядок, вы можете, наконец, перенести заказ обратно с помощью res.transpose(2,0,1) .