Интерполяция по нерегулярной сетке

Ответ 1

Попробуйте комбинацию взвешивания с обратным расстоянием и scipy.spatial.KDTree описанный в SO inverse-distance-weighted-idw-interpolation-with-python. Kd-деревья работать красиво в 2d 3d..., взвешивание с обратным расстоянием является гладким и локальным, и k = число ближайших соседей может быть изменено до скорости/точности компромисса.

Ответ 2

отличный пример обратного расстояния Роджером Вечаной и Ровирой вместе с некоторым кодом, использующим GDAL для записи в geotiff, если вы в это.

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

Копия его алгоритма с тестом:

from math import pow  
from math import sqrt  
import numpy as np  
import matplotlib.pyplot as plt  

def pointValue(x,y,power,smoothing,xv,yv,values):  
    nominator=0  
    denominator=0  
    for i in range(0,len(values)):  
        dist = sqrt((x-xv[i])*(x-xv[i])+(y-yv[i])*(y-yv[i])+smoothing*smoothing);  
        #If the point is really close to one of the data points, return the data point value to avoid singularities  
        if(dist<0.0000000001):  
            return values[i]  
        nominator=nominator+(values[i]/pow(dist,power))  
        denominator=denominator+(1/pow(dist,power))  
    #Return NODATA if the denominator is zero  
    if denominator > 0:  
        value = nominator/denominator  
    else:  
        value = -9999  
    return value  

def invDist(xv,yv,values,xsize=100,ysize=100,power=2,smoothing=0):  
    valuesGrid = np.zeros((ysize,xsize))  
    for x in range(0,xsize):  
        for y in range(0,ysize):  
            valuesGrid[y][x] = pointValue(x,y,power,smoothing,xv,yv,values)  
    return valuesGrid  


if __name__ == "__main__":  
    power=1  
    smoothing=20  

    #Creating some data, with each coodinate and the values stored in separated lists  
    xv = [10,60,40,70,10,50,20,70,30,60]  
    yv = [10,20,30,30,40,50,60,70,80,90]  
    values = [1,2,2,3,4,6,7,7,8,10]  

    #Creating the output grid (100x100, in the example)  
    ti = np.linspace(0, 100, 100)  
    XI, YI = np.meshgrid(ti, ti)  

    #Creating the interpolation function and populating the output matrix value  
    ZI = invDist(xv,yv,values,100,100,power,smoothing)  


    # Plotting the result  
    n = plt.normalize(0.0, 100.0)  
    plt.subplot(1, 1, 1)  
    plt.pcolor(XI, YI, ZI)  
    plt.scatter(xv, yv, 100, values)  
    plt.title('Inv dist interpolation - power: ' + str(power) + ' smoothing: ' + str(smoothing))  
    plt.xlim(0, 100)  
    plt.ylim(0, 100)  
    plt.colorbar()  

    plt.show() 

Ответ 3

Здесь есть несколько вариантов, которые лучше всего будут зависеть от ваших данных... Однако я не знаю, какое из готовых решений для вас

Вы говорите, что ваши входные данные взяты из трехполярных данных. Существует три основных случая, по которым эти данные могут быть структурированы.

• Отбирается из трехмерной сетки в трехполюсном пространстве, проецируется обратно на 2d LAT, данные LON.
• Отбирается из 2d-сетки в трехполюсном пространстве, проецируется в 2d LAT-данные LON.
• Неструктурированные данные в трехполюсном пространстве, проецируемые в 2d данные LAT LON

Самый простой из них - 2. Вместо интерполирования в пространстве LAT LON "просто" преобразует вашу точку обратно в исходное пространство и интерполирует там.

Другим вариантом, который работает для 1 и 2, является поиск ячеек, которые отображаются из трехполярного пространства, чтобы покрыть ваш образец. (Вы можете использовать структуру BSP или тип сетки, чтобы ускорить этот поиск). Выберите одну из ячеек и выполните интерполяцию внутри нее.

Наконец, есть куча вариантов неструктурированной интерполяции.. но они, как правило, медленные. Мой личный фаворит заключается в использовании линейной интерполяции ближайших N точек, нахождение тех N точек снова может быть выполнено с помощью грида или BSP. Другим хорошим вариантом является то, что Делоне триангулирует неструктурированные точки и интерполирует на полученную треугольную сетку.

Лично, если мой сет был случай 1, я бы использовал неструктурированную стратегию, так как я был бы обеспокоен необходимостью обрабатывать поиск через ячейки с перекрывающимися проекциями. Выбор "правильной" ячейки будет затруднен.

Ответ 4

Я предлагаю вам взглянуть на функции интерполяции GRASS (с открытым исходным кодом GIS) (http://grass.ibiblio.org/gdp/html_grass62/v.surf.bspline.html). Это не в python, но вы можете переопределить его или интерфейс с C-кодом.

Ответ 5

Я правильно понимаю, что ваши сетки данных выглядят примерно так (красный - это старые данные, синий - это новые интерполированные данные)?

alt text http://www.geekops.co.uk/photos/0000-00-02%20%28Forum%20images%29/DataSeparation.png

Это может быть немного грубый-силовой подход, но как насчет рендеринга ваших существующих данных в виде растрового изображения (opengl будет делать простую интерполяцию цветов для вас с нужными параметрами, и вы можете отображать данные как треугольники, которые должны быть довольно быстрым). Затем вы можете пробовать пиксели в местах расположения новых точек.

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

Ответ 6

Существует библиотека FORTRAN, называемая BIVAR, что очень подходит для этой проблемы. С помощью нескольких модификаций вы можете использовать его в python с помощью f2py.

Из описания:

BIVAR - библиотека FORTRAN90, которая интерполирует разбросанные двумерные данные Хироши Акимой.

BIVAR принимает набор (X, Y) точек данных, разбросанных по 2D, с соответствующими значениями данных Z и способен построить гладкую интерполяционную функцию Z (X, Y), которая согласуется с данными, и может оцениваться в других точках плоскости.