Dwuwymiarowa strukturalna interpolacja dużej tablicy z wartościami NaN lub maską

Question

Próbuję interpolować regularnie siatkowe dane windstress przy użyciu klasy Scipy's RectBivariateSpline. W niektórych punktach siatki dane wejściowe zawierają nieprawidłowe wpisy danych, które są ustawione na wartości NaN. Na początek użyłem rozwiązania do Scott's question w interpolacji dwuwymiarowej. Korzystając z moich danych, interpolacja zwraca tablicę zawierającą tylko NaN. Próbowałem również inne podejście zakładając, że moje dane są nieustrukturyzowane i przy użyciu klasy SmoothBivariateSpline. Najwyraźniej mam zbyt wiele punktów danych, aby użyć interpolacji nieustrukturyzowanej, ponieważ kształt macierzy danych wynosi (719 x 2880).

Aby zilustrować mój problem stworzyłem następujący skrypt:

from __future__ import division 
import numpy 
import pylab 

from scipy import interpolate 


# The signal and lots of noise 
M, N = 20, 30 # The shape of the data array 
y, x = numpy.mgrid[0:M+1, 0:N+1] 
signal = -10 * numpy.cos(x/50 + y/10)/(y + 1) 
noise = numpy.random.normal(size=(M+1, N+1)) 
z = signal + noise 


# Some holes in my dataset 
z[1:2, 0:2] = numpy.nan 
z[1:2, 9:11] = numpy.nan 
z[0:1, :12] = numpy.nan 
z[10:12, 17:19] = numpy.nan 


# Interpolation! 
Y, X = numpy.mgrid[0.125:M:0.5, 0.125:N:0.5] 
sp = interpolate.RectBivariateSpline(y[:, 0], x[0, :], z) 
Z = sp(Y[:, 0], X[0, :]) 

sel = ~numpy.isnan(z) 
esp = interpolate.SmoothBivariateSpline(y[sel], x[sel], z[sel], 0*z[sel]+5) 
eZ = esp(Y[:, 0], X[0, :]) 


# Comparing the results 
pylab.close('all') 
pylab.ion() 

bbox = dict(edgecolor='w', facecolor='w', alpha=0.9) 
crange = numpy.arange(-15., 16., 1.) 

fig = pylab.figure() 
ax = fig.add_subplot(1, 3, 1) 
ax.contourf(x, y, z, crange) 
ax.set_title('Original') 
ax.text(0.05, 0.98, 'a)', ha='left', va='top', transform=ax.transAxes, 
    bbox=bbox) 

bx = fig.add_subplot(1, 3, 2, sharex=ax, sharey=ax) 
bx.contourf(X, Y, Z, crange) 
bx.set_title('Spline') 
bx.text(0.05, 0.98, 'b)', ha='left', va='top', transform=bx.transAxes, 
    bbox=bbox) 

cx = fig.add_subplot(1, 3, 3, sharex=ax, sharey=ax) 
cx.contourf(X, Y, eZ, crange) 
cx.set_title('Expected') 
cx.text(0.05, 0.98, 'c)', ha='left', va='top', transform=cx.transAxes, 
    bbox=bbox) 

co daje następujące wyniki:

Rysunek przedstawia wykonaną mapę danych (a), a wyniki wykorzystujące scipy za RectBivariateSpline (b) i klasy SmoothBivariateSpline (c). Pierwsza interpolacja skutkuje tablicą zawierającą tylko NaN. Idealnie spodziewam się wyniku podobnego do drugiej interpolacji, która jest bardziej intensywna obliczeniowo. Nie potrzebuję ekstrapolacji danych poza region domeny.

Answer 1

Można użyć griddata, jedynym problemem jest to, że nie radzi sobie dobrze z krawędziami. To może być wspomagany przez na przykład odzwierciedla, w zależności od tego, jak dane są ... Oto przykład:

from __future__ import division 
import numpy 
import pylab 
from scipy import interpolate 


# The signal and lots of noise 
M, N = 20, 30 # The shape of the data array 
y, x = numpy.mgrid[0:M+1, 0:N+1] 
signal = -10 * numpy.cos(x/50 + y/10)/(y + 1) 
noise = numpy.random.normal(size=(M+1, N+1)) 
z = signal + noise 


# Some holes in my dataset 
z[1:2, 0:2] = numpy.nan 
z[1:2, 9:11] = numpy.nan 
z[0:1, :12] = numpy.nan 
z[10:12, 17:19] = numpy.nan 

zi = numpy.vstack((z[::-1,:],z)) 
zi = numpy.hstack((zi[:,::-1], zi)) 
y, x = numpy.mgrid[0:2*(M+1), 0:2*(N+1)] 
y *= 5 # anisotropic interpolation if needed. 

zi = interpolate.griddata((y[~numpy.isnan(zi)], x[~numpy.isnan(zi)]), 
       zi[~numpy.isnan(zi)], (y, x), method='cubic') 
zi = zi[:(M+1),:(N+1)][::-1,::-1] 

pylab.subplot(1,2,1) 
pylab.imshow(z, origin='lower') 
pylab.subplot(1,2,2) 
pylab.imshow(zi, origin='lower') 
pylab.show() 

Jeśli zabraknie pamięci, można podzielić dane, wzdłuż linii:

def large_griddata(data_x, vals, grid, method='nearest'): 
    x, y = data_x 
    X, Y = grid 
    try: 
     return interpolate.griddata((x,y),vals,(X,Y),method=method) 
    except MemoryError: 
     pass 

    N = (np.min(X)+np.max(X))/2. 
    M = (np.min(Y)+np.max(Y))/2. 

    masks = [(x<N) & (y<M), 
      (x<N) & (y>=M), 
      (x>=N) & (y<M), 
      (x>=N) & (y>=M)] 

    grid_mask = [(X<N) & (Y<M), 
      (X<N) & (Y>=M), 
      (X>=N) & (Y<M), 
      (X>=N) & (Y>=M)] 

    Z = np.zeros_like(X) 
    for i in range(4): 
     Z[grid_mask[i]] = large_griddata((x[masks[i]], y[masks[i]]), 
        vals[masks[i]], (X[grid_mask[i]], Y[grid_mask[i]]), method=method) 

    return Z