¿Cómo puedo realizar una interpolación bidimensional usando scipy?

Esta Q&A está pensada como una interpolación bidimensional (y multidimensional) canónica (-ish) con scipy. A menudo hay preguntas sobre la syntax básica de varios métodos de interpolación multidimensional, espero aclarar esto también.

Tengo un conjunto de puntos de datos bidimensionales dispersos, y me gustaría contourf como una buena superficie, preferiblemente usando algo como contourf o plot_surface en matplotlib.pyplot . ¿Cómo puedo interpolar mis datos bidimensionales o multidimensionales en una malla usando scipy?

He encontrado el scipy.interpolate , pero sigo recibiendo errores cuando uso interp2d o bisplrep o griddata o rbf . ¿Cuál es la syntax adecuada de estos métodos?

Descargo de responsabilidad: principalmente escribo esta publicación teniendo en cuenta las consideraciones sintácticas y el comportamiento general. No estoy familiarizado con la memoria y el aspecto de la CPU de los métodos descritos, y dirijo esta respuesta a aquellos que tienen conjuntos de datos razonablemente pequeños, de modo que la calidad de la interpolación pueda ser el aspecto principal a considerar. Soy consciente de que cuando se trabaja con conjuntos de datos muy grandes, los métodos de mejor rendimiento (a saber, griddata y Rbf ) podrían no ser factibles.

Voy a comparar tres tipos de métodos de interpolación multidimensional ( interp2d / splines, griddata y Rbf ). Los someteré a dos tipos de tareas de interpolación y dos tipos de funciones subyacentes (puntos a partir de los cuales se deben interpolar). Los ejemplos específicos demostrarán la interpolación bidimensional, pero los métodos viables son aplicables en dimensiones arbitrarias. Cada método proporciona varios tipos de interpolación; en todos los casos usaré interpolación cúbica (o algo cerca 1 ). Es importante tener en cuenta que siempre que utilice la interpolación, introducirá un sesgo en comparación con sus datos sin procesar, y los métodos específicos utilizados afectarán los artefactos con los que terminará. Siempre se consciente de esto, e interpolar responsablemente.

Las dos tareas de interpolación serán

  1. muestreo (los datos de entrada están en una cuadrícula rectangular, los datos de salida están en una cuadrícula más densa)
  2. interpolación de datos dispersos en una cuadrícula regular

Las dos funciones (sobre el dominio [x,y] in [-1,1]x[-1,1] ) serán

  1. una función suave y amigable: cos(pi*x)*sin(pi*y) ; rango en [-1, 1]
  2. una función maligna (y en particular, no continua): x*y/(x^2+y^2) con un valor de 0.5 cerca del origen; rango en [-0.5, 0.5]

Así es como se ven:

fig1: funciones de prueba

Primero demostraré cómo se comportan los tres métodos en estas cuatro pruebas, luego detallaré la syntax de las tres. Si sabe lo que debe esperar de un método, es posible que no desee perder su tiempo aprendiendo su syntax (mirándole a usted, interp2d ).

Datos de prueba

En aras de la explicación, aquí está el código con el que generé los datos de entrada. Mientras que en este caso específico obviamente estoy al tanto de la función que subyace a los datos, solo usaré esto para generar información para los métodos de interpolación. Uso numpy por conveniencia (y principalmente para generar los datos), pero bastaría con scipy solo.

 import numpy as np import scipy.interpolate as interp # auxiliary function for mesh generation def gimme_mesh(n): minval = -1 maxval = 1 # produce an asymmetric shape in order to catch issues with transpositions return np.meshgrid(np.linspace(minval,maxval,n), np.linspace(minval,maxval,n+1)) # set up underlying test functions, vectorized def fun_smooth(x, y): return np.cos(np.pi*x)*np.sin(np.pi*y) def fun_evil(x, y): # watch out for singular origin; function has no unique limit there return np.where(x**2+y**2>1e-10, x*y/(x**2+y**2), 0.5) # sparse input mesh, 6x7 in shape N_sparse = 6 x_sparse,y_sparse = gimme_mesh(N_sparse) z_sparse_smooth = fun_smooth(x_sparse, y_sparse) z_sparse_evil = fun_evil(x_sparse, y_sparse) # scattered input points, 10^2 altogether (shape (100,)) N_scattered = 10 x_scattered,y_scattered = np.random.rand(2,N_scattered**2)*2 - 1 z_scattered_smooth = fun_smooth(x_scattered, y_scattered) z_scattered_evil = fun_evil(x_scattered, y_scattered) # dense output mesh, 20x21 in shape N_dense = 20 x_dense,y_dense = gimme_mesh(N_dense) 

Función suave y upmpling

Vamos a empezar con la tarea más fácil. A continuación, se muestra cómo funciona el muestreo desde una malla de forma [6,7] a una de [20,21] para la función de prueba suave:

fig2: muestreo suave

Aunque esta es una tarea simple, ya hay diferencias sutiles entre las salidas. A primera vista, las tres salidas son razonables. Hay dos características a tener en cuenta, basadas en nuestro conocimiento previo de la función subyacente: el caso intermedio de griddata el que más distorsiona los datos. Tenga en cuenta el límite y==-1 del trazado (más cercano a la etiqueta x ): la función debe ser estrictamente cero (ya que y==-1 es una línea nodal para la función suave), sin embargo, este no es el caso de los griddata de griddata . También tenga en cuenta el límite x==-1 de los gráficos (detrás, a la izquierda): la función subyacente tiene un máximo local (que implica un gradiente cero cerca del límite) en [-1, -0.5] , pero la salida de datos de griddata muestra claramente gradiente distinto de cero en esta región. El efecto es sutil, pero es un sesgo, no obstante. (La fidelidad de Rbf es aún mejor con la opción predeterminada de funciones radiales, doblada como multiquadratic ).

Mal funcionamiento y upmpling.

Una tarea un poco más difícil es realizar muestreos en nuestra función maligna:

fig3: malversación

Se están empezando a mostrar diferencias claras entre los tres métodos. Al observar las plots de superficie, hay claros extremos espurios que aparecen en la salida de interp2d (observe las dos jorobas en el lado derecho de la superficie trazada). Si bien griddata y Rbf parecen producir resultados similares a primera vista, este último parece producir un mínimo más profundo cerca de [0.4, -0.4] que está ausente de la función subyacente.

Sin embargo, hay un aspecto crucial en el que Rbf es muy superior: respeta la simetría de la función subyacente (que, por supuesto, también es posible gracias a la simetría de la malla de la muestra). La salida de griddata rompe la simetría de los puntos de muestra, que ya es débilmente visible en el caso liso.

Función suave y datos dispersos.

La mayoría de las veces, uno quiere realizar la interpolación en datos dispersos. Por esta razón espero que estas pruebas sean más importantes. Como se muestra arriba, los puntos de muestra se eligieron pseudo-uniformemente en el dominio de interés. En escenarios realistas, es posible que tenga ruido adicional con cada medición, y debe considerar si tiene sentido interpolar sus datos en bruto para empezar.

Salida para la función lisa:

fig4: suave interpolación dispersada

Ahora ya hay un pequeño espectáculo de terror. interp2d la salida de interp2d a entre [-1, 1] exclusivamente para el trazado, a fin de preservar al menos una cantidad mínima de información. Está claro que, si bien algunas de las formas subyacentes están presentes, existen enormes regiones ruidosas donde el método se descompone completamente. El segundo caso de griddata reproduce la forma bastante bien, pero tenga en cuenta las regiones blancas en el borde del gráfico de contorno. Esto se debe al hecho de que los datos de griddata solo funcionan dentro del casco convexo de los puntos de datos de entrada (en otras palabras, no se realiza ninguna extrapolación ). Mantuve el valor predeterminado de NaN para los puntos de salida que se encuentran fuera del casco convexo. 2 Teniendo en cuenta estas características, Rbf parece tener el mejor rendimiento.

Función maligna y datos dispersos.

Y el momento que todos hemos estado esperando:

fig5: mal interpolación dispersada

No es una gran sorpresa que interp2d . De hecho, durante la llamada a interp2d debe esperar que algunos RuntimeWarning s RuntimeWarning quejen de la imposibilidad de construir la spline. En cuanto a los otros dos métodos, Rbf parece producir la mejor salida, incluso cerca de los límites del dominio donde se extrapola el resultado.


Así que permítanme decir algunas palabras acerca de los tres métodos, en orden decreciente de preferencia (de modo que lo peor sea lo menos probable que alguien lea).

scipy.interpolate.Rbf

La clase Rbf significa “funciones de base radial”. Para ser honesto, nunca he considerado este enfoque hasta que empecé a investigar para esta publicación, pero estoy bastante seguro de que los usaré en el futuro.

Al igual que los métodos basados ​​en spline (ver más adelante), el uso viene en dos pasos: el primero crea una instancia de clase Rbf se puede llamar en base a los datos de entrada, y luego llama a este objeto para obtener una malla de salida determinada para obtener el resultado interpolado. Ejemplo de la prueba de muestreo suave:

 import scipy.interpolate as interp zfun_smooth_rbf = interp.Rbf(x_sparse, y_sparse, z_sparse_smooth, function='cubic', smooth=0) # default smooth=0 for interpolation z_dense_smooth_rbf = zfun_smooth_rbf(x_dense, y_dense) # not really a function, but a callable class instance 

Tenga en cuenta que tanto los puntos de entrada como los de salida eran matrices 2d en este caso, y la salida z_dense_smooth_rbf tiene la misma forma que x_dense y y_dense sin ningún esfuerzo. También tenga en cuenta que Rbf admite dimensiones arbitrarias para la interpolación.

Por lo tanto, scipy.interpolate.Rbf

  • produce una salida de buen comportamiento incluso para datos de entrada locos
  • Soporta interpolación en dimensiones superiores.
  • extrapola fuera del casco convexo de los puntos de entrada (por supuesto, la extrapolación siempre es un juego de azar, y por lo general no debe confiar en ello)
  • crea un interpolador como primer paso, por lo que evaluarlo en varios puntos de salida es menos esfuerzo adicional
  • puede tener puntos de salida de forma arbitraria (en lugar de estar restringido a mallas rectangulares, ver más adelante)
  • propensos a preservar la simetría de los datos de entrada
  • soporta múltiples tipos de funciones radiales para la función de palabras clave: multiquadric , inverse , gaussian , linear , cubic , thin_plate , thin_plate y arbitrario definido por el usuario

scipy.interpolate.griddata

Mi anterior favorito, griddata , es un caballo de batalla general para la interpolación en dimensiones arbitrarias. No realiza extrapolaciones más allá de establecer un único valor preestablecido para puntos fuera del casco convexo de los puntos nodales, pero dado que la extrapolación es algo muy voluble y peligroso, esto no es necesariamente una estafa. Ejemplo de uso:

 z_dense_smooth_griddata = interp.griddata(np.array([x_sparse.ravel(),y_sparse.ravel()]).T, z_sparse_smooth.ravel(), (x_dense,y_dense), method='cubic') # default method is linear 

Tenga en cuenta la syntax ligeramente kludgy. Los puntos de entrada deben especificarse en una matriz de forma [N, D] en dimensiones D Para esto, primero tenemos que aplanar nuestras matrices de coordenadas 2d (utilizando ravel ), luego concatenar las matrices y transponer el resultado. Hay varias formas de hacer esto, pero todas parecen ser voluminosas. Los datos de entrada z también tienen que ser aplanados. Tenemos un poco más de libertad cuando se trata de los puntos de salida: por alguna razón, estos también pueden especificarse como una tupla de matrices multidimensionales. Tenga en cuenta que la help de griddata es engañosa, ya que sugiere que lo mismo es cierto para los puntos de entrada (al menos para la versión 0.17.0):

 griddata(points, values, xi, method='linear', fill_value=nan, rescale=False) Interpolate unstructured D-dimensional data. Parameters ---------- points : ndarray of floats, shape (n, D) Data point coordinates. Can either be an array of shape (n, D), or a tuple of `ndim` arrays. values : ndarray of float or complex, shape (n,) Data values. xi : ndarray of float, shape (M, D) Points at which to interpolate data. 

En pocas palabras, scipy.interpolate.griddata

  • produce una salida de buen comportamiento incluso para datos de entrada locos
  • Soporta interpolación en dimensiones superiores.
  • no realiza la extrapolación, se puede establecer un solo valor para la salida fuera del casco convexo de los puntos de entrada (ver fill_value )
  • calcula los valores interpolados en una sola llamada, por lo que el sondeo de múltiples conjuntos de puntos de salida comienza desde cero
  • Puede tener puntos de salida de forma arbitraria.
  • soporta interpolación lineal y vecina más cercana en dimensiones arbitrarias, cúbica en 1d y 2d. La interpolación lineal más cercana y cercana utiliza el Interpolador NearestNDInterpolator más cercano y el Interpolador NearestNDInterpolator LinearNDInterpolator debajo de la cubierta, respectivamente. La interpolación cúbica 1d usa una spline, la interpolación cúbica 2d usa el Interpolador CloughTocher2DInterpolator para construir un interpolador cúbico en partes continuamente diferenciable.
  • podría violar la simetría de los datos de entrada

scipy.interpolate.interp2d / scipy.interpolate.bisplrep

La única razón por la que interp2d y sus familiares es que tiene un nombre engañoso y es probable que la gente intente usarlo. Alerta de spoiler: no lo use (a partir de la versión 0.17.0 de scipy). Ya es más especial que los temas anteriores, ya que se usa específicamente para la interpolación bidimensional, pero sospecho que este es el caso más común para la interpolación multivariable.

En lo que respecta a la syntax, interp2d es similar a Rbf que primero necesita construir una instancia de interpolación, que puede llamarse para proporcionar los valores interpolados reales. Sin embargo, hay un problema: los puntos de salida deben ubicarse en una malla rectangular, por lo que las entradas que entran en la llamada al interpolador tienen que ser vectores 1d que abarquen la cuadrícula de salida, como desde numpy.meshgrid :

 # reminder: x_sparse and y_sparse are of shape [6, 7] from numpy.meshgrid zfun_smooth_interp2d = interp.interp2d(x_sparse, y_sparse, z_sparse_smooth, kind='cubic') # default kind is 'linear' # reminder: x_dense and y_dense are of shape [20, 21] from numpy.meshgrid xvec = x_dense[0,:] # 1d array of unique x values, 20 elements yvec = y_dense[:,0] # 1d array of unique y values, 21 elements z_dense_smooth_interp2d = zfun_smooth_interp2d(xvec,yvec) # output is [20, 21]-shaped array 

Uno de los errores más comunes cuando se usa interp2d es poner sus mallas 2d completas en la llamada de interpolación, lo que lleva a un consumo explosivo de memoria, y con suerte a un MemoryError apresurado.

Ahora, el mayor problema con interp2d es que a menudo no funciona. Para entender esto, tenemos que mirar debajo del capó. Resulta que interp2d es una envoltura para las funciones de nivel inferior bisplrep + bisplev , que a su vez son envolturas para rutinas FITPACK (escritas en Fortran). La llamada equivalente al ejemplo anterior sería

 kind = 'cubic' if kind=='linear': kx=ky=1 elif kind=='cubic': kx=ky=3 elif kind=='quintic': kx=ky=5 # bisplrep constructs a spline representation, bisplev evaluates the spline at given points bisp_smooth = interp.bisplrep(x_sparse.ravel(),y_sparse.ravel(),z_sparse_smooth.ravel(),kx=kx,ky=ky,s=0) z_dense_smooth_bisplrep = interp.bisplev(xvec,yvec,bisp_smooth).T # note the transpose 

Ahora, esto es lo que tiene que ver con interp2d : (en scipy versión 0.17.0) hay un buen comentario en interpolate/interpolate.py para interp2d :

 if not rectangular_grid: # TODO: surfit is really not meant for interpolation! self.tck = fitpack.bisplrep(x, y, z, kx=kx, ky=ky, s=0.0) 

y de hecho, en interpolate/fitpack.py , en bisplrep hay alguna configuración y, en última instancia,

 tx, ty, c, o = _fitpack._surfit(x, y, z, w, xb, xe, yb, ye, kx, ky, task, s, eps, tx, ty, nxest, nyest, wrk, lwrk1, lwrk2) 

Y eso es. Las rutinas que subyacen a interp2d no están realmente destinadas a realizar la interpolación. Pueden ser suficientes para datos suficientemente bien comportados, pero en circunstancias realistas probablemente querrá usar otra cosa.

Solo para concluir, interpolate.interp2d

  • puede llevar a artefactos incluso con datos bien temperados
  • es específicamente para problemas bivariados (aunque hay una interpn limitada para los puntos de entrada definidos en una cuadrícula)
  • realiza extrapolación
  • crea un interpolador como primer paso, por lo que evaluarlo en varios puntos de salida es menos esfuerzo adicional
  • solo puede producir una salida sobre una cuadrícula rectangular, para una salida dispersa tendría que llamar al interpolador en un bucle
  • Soporta interpolación lineal, cúbica y quíntica.
  • podría violar la simetría de los datos de entrada

1 Estoy bastante seguro de que las funciones de base cubic y linear de Rbf no se corresponden exactamente con los otros interpoladores del mismo nombre.
2 Estos NaN también son la razón por la cual el diagtwig de superficie parece tan extraño: matplotlib históricamente tiene dificultades para trazar objetos complejos en 3D con información de profundidad adecuada. Los valores de NaN en los datos confunden al renderizador, por lo que las partes de la superficie que deberían estar en la parte posterior se trazan en la parte frontal. Este es un problema con la visualización, y no con la interpolación.