Interpolación bivariada estructurada de una matriz grande con valores o máscara de NaN

Estoy tratando de interpolar los datos del parabrisas con una rejilla RectBivariateSpline de Scipy’s. En algunos puntos de la cuadrícula, los datos de entrada contienen entradas de datos no válidas, que se establecen en valores NaN. Para empezar, utilicé la solución a la pregunta de Scott sobre la interpolación bidimensional. Usando mis datos, la interpolación devuelve una matriz que contiene solo NaNs. También he intentado un enfoque diferente asumiendo que mis datos no están estructurados y que usan la clase SmoothBivariateSpline . Aparentemente tengo demasiados puntos de datos para usar la interpolación no estructurada, ya que la forma de la matriz de datos es (719 x 2880).

Para ilustrar mi problema he creado el siguiente 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) 

Lo que da los siguientes resultados: Rejillas bivariadas. (a) Los datos originales construidos, (b) las clases RectBivariateSpline de Scipy y (c) las clases SmoothBivariateSpline.

La figura muestra un mapa de datos construido (a) y los resultados utilizando las clases RectBivariateSpline (b) y SmoothBivariateSpline (c) de Scipy. La primera interpolación da como resultado una matriz con solo NaN. Idealmente, habría esperado un resultado similar a la segunda interpolación, que es más intensivo computacionalmente. No necesariamente necesito extrapolación de datos fuera de la región del dominio.

Podría usar griddata , el único problema es que no maneja bien los bordes. Esto podría ser de ayuda, por ejemplo, reflexionando, dependiendo de cómo estén sus datos … Aquí hay un ejemplo:

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

Si se queda sin memoria, puede dividir sus datos, en la línea de:

 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=M), (x>=N) & (y=N) & (y>=M)] grid_mask = [(X=M), (X>=N) & (Y=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