Как интерполировать / экстраполировать в частично пустой регулярной сетке? - PullRequest
Новая
       7

Как интерполировать / экстраполировать в частично пустой регулярной сетке?

0 голосов
/ 02 мая 2018

Я хотел бы создать функцию python для линейной интерполяции в частично пустой сетке и получения ближайшей экстраполяции за пределами границ.

Допустим, у меня есть следующие данные, хранящиеся в пандах DataFrame: 

In [1]: import numpy as np
In [2]: import pandas as pd

In [3]: x = [0,1,2,3,4]
In [4]: y = [0.5,1.5,2.5,3.5,4.5,5.5]
In [5]: z = np.array([[np.nan,np.nan,1.5,2.0,5.5,3.5],[np.nan,1.0,4.0,2.5,4.5,3.0],[2.0,0.5,6.0,1.5,3.5,np.nan],[np.nan,1.5,4.0,2.0,np.nan,np.nan],[np.nan,np.nan,2.0,np.nan,np.nan,np.nan]])
In [6]: df = pd.DataFrame(t,index=x,columns=y)
In [7]: df
Out[7]:
    0.5  1.5  2.5  3.5  4.5  5.5
 0  NaN  NaN  1.5  2.0  5.5  3.5
 1  NaN  1.0  4.0  2.5  4.5  3.0
 2  2.0  0.5  6.0  1.5  3.5  NaN
 3  NaN  1.5  4.0  2.0  NaN  NaN
 4  NaN  NaN  2.0  NaN  NaN  NaN

Я хотел бы получить функцию myInterp, которая возвращает линейную интерполяцию в пределах границ данных (т. Е. Не значений NaN) и получать ближайшую экстраполяцию за пределами границ (т. Е. NaN или без значений), например: 

In [1]: myInterp([1.5,2.5]) #linear interpolation
Out[1]: 5.0

In [2]: myInterp([1.5,4.0]) #bi-linear interpolation
Out[2]: 3.0

In [3]: myInterp([0.0,2.0]) #nearest extrapolation (inside grid)
Out[3]: 1.5

In [4]: myInterp([5.0,2.5]) #nearest extrapolation (outside grid)
Out[4]: 2.0

Я пробовал многие комбинации scipy.interpolate безуспешно, у кого-нибудь есть предложения, как это сделать?

1 Ответ

0 голосов
/ 02 мая 2018

Да, к сожалению, Сципи не имеет дела с Нансом

Из документов: 

Note that calling interp2d with NaNs present in input values results in undefined behaviour.

Даже маскировка nans в np.masked_array не удалась.

Поэтому мой совет - удалить все записи nan из z, воспользовавшись возможностью, чтобы дать sp.interp2d полный список координат x и y только для действительных данных и оставить z также 1D: 

X=[];Y=[];Z=[]                     # initialize new 1-D-lists for interp2
for i, xi in enumerate(x):         # iterate through x
    for k, yk in enumerate(y):     # iterate through y
        if not np.isnan(z[i, k]):  # check if z-value is valid...
            X.append(xi)           # ...and if so, append coordinates and value to prepared lists
            Y.append(yk)
            Z.append(z[i, k])

Таким образом, по крайней мере, sp.interp2d работает и дает результат: 

ip = sp.interpolate.interp2d(X,Y,Z)

Однако значения в результате вас не порадуют: 

In: ip(x,y)
Out: 
array([[ 18.03583061,  -0.44933642,   0.83333333,  -1.        , -1.46105542],
       [  9.76791531,   1.3014037 ,   2.83333333,   1.5       ,  0.26947229],
       [  1.5       ,   3.05214381,   4.83333333,   4.        ,   2.        ],
       [  2.        ,   3.78378051,   1.5       ,   2.        ,   0.8364618 ],
       [  5.5       ,   3.57039277,   3.5       ,  -0.83019815,  -0.7967441 ],
       [  3.5       ,   3.29227922,  17.29607177,   0.        ,   0.        ]])

по сравнению с входными данными: 

In:z
Out: 
array([[ nan,  nan,  1.5,  2. ,  5.5,  3.5],
       [ nan,  1. ,  4. ,  2.5,  4.5,  3. ],
       [ 2. ,  0.5,  6. ,  1.5,  3.5,  nan],
       [ nan,  1.5,  4. ,  2. ,  nan,  nan],
       [ nan,  nan,  2. ,  nan,  nan,  nan]])

Но ИМХО это потому, что изменения градиента в ваших данных слишком велики. Еще больше относительно небольшого количества выборок данных.

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

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

In:ip(x,y)
Out: 
array([[ 3.        ,  3.        ,  3.        ,  3.        ,  0.        ],
       [ 3.        ,  3.        ,  3.        ,  3.        ,  1.94701008],
       [ 3.        ,  3.        ,  3.        ,  3.        ,  3.        ],
       [ 3.        ,  3.        ,  3.        ,  3.        ,  1.54973345],
       [ 3.        ,  3.        ,  3.        ,  3.        ,  0.37706713],
       [ 3.        ,  3.        ,  2.32108317,  0.75435203,  0.        ]])

результат тривиального тестового ввода 

In:z
Out: 
array([[ nan,  nan,   3.,   3.,   3.,   3.],
       [ nan,   3.,   3.,  nan,   3.,   3.],
       [  3.,   3.,   3.,   3.,   3.,  nan],
       [ nan,   3.,   3.,   3.,  nan,  nan],
       [ nan,  nan,   3.,  nan,  nan,  nan]])

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

Но сюрприз: кубическая версия работает намного лучше: 

In:ip = sp.interpolate.interp2d(X,Y,Z, kind='cubic')

In:ip(x,y)
Out: 
array([[ 3.        ,  3.        ,  3.        ,  3.02397028,  3.0958811 ],
       [ 3.        ,  3.        ,  3.        ,  3.        ,  3.        ],
       [ 3.        ,  3.        ,  3.        ,  3.        ,  3.        ],
       [ 3.        ,  3.        ,  3.        ,  3.        ,  3.        ],
       [ 3.        ,  3.        ,  3.        ,  2.97602972,  2.9041189 ],
       [ 3.        ,  3.        ,  3.        ,  2.9041189 ,  2.61647559]])

In:z
Out: 
array([[ nan,  nan,   3.,   3.,   3.,   3.],
       [ nan,   3.,   3.,  nan,   3.,   3.],
       [  3.,   3.,   3.,   3.,   3.,  nan],
       [ nan,   3.,   3.,   3.,  nan,  nan],
       [ nan,  nan,   3.,  nan,  nan,  nan]])
...