signal lagrange interpolate interpld from extrapolate python scipy interpolation

python - lagrange - ¿Cómo puedo realizar una interpolación bidimensional con scipy?



signal interpolation python (1)

Este Q&A está pensado como un canónico (-ish) con respecto a la interpolación bidimensional (y multidimensional) usando scipy. A menudo hay preguntas sobre la sintaxis básica de varios métodos de interpolación multidimensional, espero aclararlos también.

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

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

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

Voy a comparar tres tipos de métodos de interpolación multidimensionales ( interp2d / splines, griddata y Rbf ). Los someteré a dos tipos de tareas de interpolación y dos tipos de funciones subyacentes (puntos desde 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 parecido a 1 ). Es importante tener en cuenta que cada vez que utiliza la interpolación introduce sesgo en comparación con sus datos sin procesar, y los métodos específicos utilizados afectan los artefactos con los que terminará. Tenga siempre en cuenta esto e interpole de manera responsable.

Las dos tareas de interpolación serán

  1. muestreo ascendente (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 malvada (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:

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

Datos de prueba

En aras de lo explícito, aquí está el código con el que generé los datos de entrada. Si bien en este caso específico, obviamente, conozco la función subyacente a los datos, solo la usaré para generar información para los métodos de interpolación. Utilizo numpy por conveniencia (y principalmente para generar los datos), pero solo scipy sería suficiente.

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 muestreo superior

Comencemos con la tarea más fácil. Así es como un muestreo ascendente de una malla de forma [6,7] a una de [20,21] funciona para la función de prueba suave:

Aunque esta es una tarea simple, ya existen diferencias sutiles entre los resultados. 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 medio de griddata distorsiona más los datos. Tenga en cuenta el límite y==-1 del gráfico (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 para los griddata de griddata . También tenga en cuenta el límite x==-1 de las parcelas (detrás, a la izquierda): la función subyacente tiene un máximo local (lo que implica un gradiente cero cerca del límite) en [-1, -0.5] , pero la salida de datos de griddata se muestra claramente gradiente distinto de cero en esta región. El efecto es sutil, pero no obstante es un sesgo. (La fidelidad de Rbf es aún mejor con la elección predeterminada de funciones radiales, denominada multiquadratic ).

Mala función y muestreo ascendente

Una tarea un poco más difícil es realizar un muestreo superior de nuestra función malvada:

Se empiezan a mostrar claras diferencias entre los tres métodos. Al observar los gráficos de superficie, aparecen claros extremos espurios en la salida de interp2d (tenga en cuenta las dos jorobas en el lado derecho de la superficie trazada). Mientras que 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 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.

En la mayoría de los casos, se desea realizar una 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 fueron elegidos 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 sin procesar para comenzar.

Salida para la función suave:

Ahora ya hay un espectáculo de terror. interp2d la salida de interp2d 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 parte de la forma subyacente está presente, hay grandes regiones ruidosas donde el método se descompone por completo. El segundo caso de griddata reproduce la forma bastante bien, pero tenga en cuenta las regiones blancas en el borde de la gráfica de contorno. Esto se debe al hecho de que griddata solo funciona dentro del casco convexo de los puntos de datos de entrada (en otras palabras, no realiza ninguna extrapolación ). Mantuve el valor predeterminado de NaN para los puntos de salida que se encuentran fuera del casco convexo. 2 Considerando estas características, Rbf parece funcionar mejor.

Mala función y datos dispersos

Y el momento que todos hemos estado esperando:

No es una gran sorpresa que interp2d . De hecho, durante la llamada a interp2d , debe esperar algunos RuntimeWarning amistosos quejándose 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.

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

scipy.interpolate.Rbf

La clase Rbf significa "funciones de base radial". Para ser sincero, nunca he considerado este enfoque hasta que comencé 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: primero crea una instancia de clase Rbf invocable basada en los datos de entrada, y luego llama a este objeto para una malla de salida dada 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 los puntos de entrada y salida eran matrices 2d en este caso, y la salida z_dense_smooth_rbf tiene la misma forma que x_dense e y_dense sin ningún esfuerzo. También tenga en cuenta que Rbf admite dimensiones arbitrarias para la interpolación.

Entonces, scipy.interpolate.Rbf

  • produce una salida con buen comportamiento incluso para datos de entrada locos
  • admite interpolación en dimensiones superiores
  • extrapola fuera del casco convexo de los puntos de entrada (por supuesto, la extrapolación siempre es una apuesta y, por lo general, no debe confiar en ella)
  • 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 limitado a mallas rectangulares, ver más adelante)
  • propenso a preservar la simetría de los datos de entrada
  • admite múltiples tipos de funciones radiales para la función de palabras clave: multiquadric , inverse , gaussian , linear , cubic , quintic , thin_plate y arbitraria definida por el usuario

scipy.interpolate.griddata

Mi antiguo favorito, griddata , es un caballo de batalla general para la interpolación en dimensiones arbitrarias. No realiza extrapolación 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 sintaxis ligeramente kludgy. Los puntos de entrada deben especificarse en una matriz de formas [N, D] en dimensiones D Para esto, primero tenemos que aplanar nuestras matrices de coordenadas 2D (usando 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 deben ser aplanados. Tenemos un poco más de libertad en lo que respecta a los puntos de salida: por alguna razón, estos también se pueden especificar 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 con buen comportamiento incluso para datos de entrada locos
  • admite interpolación en dimensiones superiores
  • no realiza extrapolación, se puede establecer un único valor para la salida fuera del casco convexo de los puntos de entrada (consulte el valor de fill_value )
  • calcula los valores interpolados en una sola llamada, por lo que sondear múltiples conjuntos de puntos de salida comienza desde cero
  • puede tener puntos de salida de forma arbitraria
  • admite la interpolación lineal más cercana y vecina en dimensiones arbitrarias, cúbicas en 1d y 2d. La interpolación más cercana y la interpolación lineal utilizan NearestNDInterpolator y LinearNDInterpolator debajo del capó, respectivamente. La interpolación cúbica 1d usa una spline, la interpolación cúbica 2d usa CloughTocher2DInterpolator para construir un interpolador cúbico por 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 estoy discutiendo sobre interp2d y sus familiares es porque tiene un nombre engañoso y es probable que las personas intenten usarlo. Alerta de spoiler: no lo use (a partir de la versión scipy 0.17.0). 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, con mucho, el caso más común para la interpolación multivariante.

En cuanto a la sintaxis, interp2d es similar a Rbf en que primero necesita construir una instancia de interpolación, a la que se pueda llamar 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 deben ser vectores 1d que abarcan la cuadrícula de salida, como si numpy.meshgrid de 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 al usar interp2d es poner sus mallas 2D completas en la llamada de interpolación, lo que conduce 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 un contenedor para las funciones de nivel inferior bisplrep + bisplev , que a su vez son contenedores para las 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 pasa con interp2d : (en la versión scipy 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 finalmente

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 subyacentes a interp2d no están realmente destinadas a realizar interpolación. Pueden ser suficientes para obtener datos suficientemente bien comportados, pero en circunstancias realistas probablemente querrá usar otra cosa.

Solo para concluir, interpolate.interp2d

  • puede conducir a artefactos incluso con datos moderados
  • es específicamente para problemas bivariados (aunque existe 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 en una cuadrícula rectangular, para una salida dispersa debería llamar al interpolador en un bucle
  • admite 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 el tipo cubic y linear de funciones Rbf 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 diagrama de superficie parece tan extraño: matplotlib históricamente tiene dificultades para trazar objetos 3D complejos 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 el reverso se trazan para estar en el frente. Este es un problema con la visualización, y no con la interpolación.