Interpolazione strutturata bivariata di array di grandi dimensioni con valori NaN o maschera

Question

Sto tentando di interpolare regolarmente dati grigliati windstress utilizzando la classe RectBivariateSpline di Scipy. In alcuni punti della griglia, i dati di input contengono voci di dati non valide, che sono impostate su valori NaN. Per cominciare, ho usato la soluzione su Scott's question sull'interpolazione bidimensionale. Usando i miei dati, l'interpolazione restituisce un array contenente solo NaN. Ho anche provato un approccio diverso ipotizzando che i miei dati non siano strutturati e che utilizzi la classe SmoothBivariateSpline. Apparentemente ho troppi punti dati per usare l'interpolazione non strutturata, dato che la forma dell'array di dati è (719 x 2880).

Per illustrare il mio problema ho creato il seguente script:

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) 

che dà i seguenti risultati:

La figura mostra una mappa di dati costruito (a) ed i risultati utilizzando RectBivariateSpline di SciPy (b) e SmoothBivariateSpline (c) classi. La prima interpolazione produce un array con solo NaN. Idealmente mi sarei aspettato un risultato simile alla seconda interpolazione, che è più intensivo dal punto di vista computazionale. Non ho necessariamente bisogno di estrapolazione dei dati al di fuori della regione del dominio.

Answer 1

È possibile utilizzare griddata, l'unico problema è che non gestisce bene i bordi. Questo potrebbe essere aiutato da riflettere per esempio, a seconda di come i dati sono al ... Ecco un esempio:

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() 

Se si esaurisce la memoria, è possibile suddividere i dati, sulla falsariga di:

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