Интерполяция по неправильной сетке
Так, у меня есть три массива numpy, которые хранят широту, долготу и некоторое значение свойства на сетке - то есть, я имею LAT (y, x), LON (y, x), и, говорю температуру T (y, x), для некоторых пределов X и Y. Сетка не обязательно обычна - на самом деле, это является трехполюсным.
Я затем хочу интерполировать эти свойство (температурные) значения на набор различных точек lat/lon (сохраненный как lat1 (t), lon1 (t), приблизительно для 10 000 т...), которые не падают на фактические узлы решетки. Я попробовал matplotlib.mlab.griddata, но это берет слишком долго (он действительно не разработан для того, что я делаю, в конце концов). Я также попробовал scipy.interpolate.interp2d, но я получаю MemoryError (мои сетки о 400x400).
Есть ли какой-либо вид пленки, предпочтительно быстрый способ сделать это? Я не могу не думать, что ответ является чем-то очевидным... Спасибо!!
4 ответа
Попробуйте использовать комбинацию взвешивания с обратным расстоянием и scipy.spatial.KDTree , описанную в SO inverse-distance-weighted-idw-interpolation-with-python. Kd-trees отлично работают в 2d 3d ..., взвешивание по обратному расстоянию гладкое и локальное, а k=количество ближайших соседей можно варьировать, чтобы найти компромисс между скоростью и точностью.
Здесь есть несколько вариантов, и какой из них лучше всего будет зависеть от ваших данных ... Однако я не знаю готового решения для вас
Вы говорите, что ваши входные данные взяты из трехполюсных данных. Есть три основных случая того, как эти данные могут быть структурированы.
- Выборка из трехмерной сетки в трехполюсном пространстве, спроецированная обратно на данные 2d LAT, LON.
- Отобраны из 2-мерной сетки в трехполярном пространстве, спроецированы в 2-мерные данные LAT LON.
- Неструктурированные данные в трехполюсном пространстве проецируются в данные 2d LAT LON
Самый простой из них - 2. Вместо интерполяции в пространстве LAT LON «просто» преобразуйте вашу точку обратно в исходное пространство и интерполируйте там.
Другой вариант, который работает для 1 и 2, - это поиск ячеек, отображаемых из трехполюсного пространства для покрытия вашей точки выборки. (Вы можете использовать BSP или структуру типа сетки, чтобы ускорить этот поиск) Выберите одну из ячеек и выполните интерполяцию внутри нее.
Наконец, есть куча вариантов неструктурированной интерполяции ... но они, как правило, медленные. Мой личный фаворит - использовать линейную интерполяцию ближайших N точек, поиск этих N точек снова может быть выполнен с использованием координатной сетки или BSP. Еще один хороший вариант - триангулировать неструктурированные точки Делоне и интерполировать полученную треугольную сетку.
Лично, если бы моя сетка была случаем 1, я бы использовал неструктурированную стратегию, так как я бы беспокоился о необходимости обрабатывать поиск по ячейкам с перекрывающимися проекциями. Выбрать «правильную» ячейку будет сложно.
Я предлагаю вам взглянуть на функции интерполяции GRASS (ГИС-пакет с открытым исходным кодом) (http://grass.ibiblio.org/gdp/html_grass62/v.surf.bspline.html). Это не в python, но вы можете переделать его или сопрячь с кодом на C.
Правильно ли я понимаю, что ваши сетки данных выглядят примерно так (красный - старые данные, синий - новые интерполированные данные)?
alt text http://www.geekops.co.uk/photos/0000-00-02%20%28Forum%20images%29/DataSeparation.png
Возможно, это немного грубый подход, но как насчет рендеринга существующих данных в виде растрового изображения (opengl сделает для вас простую интерполяцию цветов при правильных настройках, и вы сможете рендерить данные в виде треугольников, что должно быть довольно быстро). Затем можно сделать выборку пикселей в местах расположения новых точек.
В качестве альтернативы можно отсортировать первый набор точек пространственно, а затем найти ближайшие старые точки, окружающие новую точку, и интерполировать на основе расстояний до этих точек.