Esta sesión de preguntas y respuestas está pensada como un canónico (-ish) relativo a la interpolación bidimensional (y multidimensional) utilizando scipy. A menudo hay preguntas sobre la sintaxis básica de varios métodos de interpolación multidimensional, espero aclararlas también.

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

Encontré el scipy.interpolatesubpaquete, pero sigo recibiendo errores al usar interp2do bisplrepo griddatao rbf. ¿Cuál es la sintaxis adecuada de estos métodos?

Descargo de responsabilidad: principalmente estoy escribiendo esta publicación teniendo en cuenta consideraciones sintácticas y el comportamiento general. No estoy familiarizado con el aspecto de memoria y 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 puede ser el aspecto principal a considerar. Soy consciente de que cuando se trabaja con conjuntos de datos muy grandes, es posible que los métodos de mejor rendimiento (a saber, griddatay Rbf) no sean factibles.

Voy a comparar tres tipos de métodos de interpolación multidimensionales ( interp2d/ splines griddatay Rbf). Los someteré a dos tipos de tareas de interpolación y dos tipos de funciones subyacentes (puntos a partir de los cuales se van a 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 cercano a ¹ ). Es importante tener en cuenta que cada vez que usa la interpolación, introduce sesgos en comparación con sus datos sin procesar, y los métodos específicos utilizados afectan los artefactos con los que terminará. Sea siempre consciente de 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 de suavizar 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:

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

Datos de prueba

En aras de la claridad, aquí está el código con el que generé los datos de entrada. Si bien en este caso específico obviamente soy consciente de la función subyacente a los datos, solo usaré esto para generar entradas para los métodos de interpolación. Uso numpy por conveniencia (y principalmente para generar los datos), pero scipy solo sería suficiente también.

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 sencilla. Así es como funciona un muestreo superior de una malla de forma [6,7]a una de [20,21]las funciones de prueba suave:

Aunque se trata de una tarea sencilla, ya existen diferencias sutiles entre los resultados. A primera vista, los tres resultados son razonables. Hay dos características a tener en cuenta, basadas en nuestro conocimiento previo de la función subyacente: el caso intermedio de griddatadistorsiona más los datos. Tenga en cuenta el y==-1límite del gráfico (más cercano a la xetiqueta): la función debe ser estrictamente cero (ya que y==-1es una línea nodal para la función suave), pero este no es el caso de griddata. También tenga en cuenta el x==-1límite de los gráficos (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], sin embargo, la griddatasalida muestra un gradiente claramente distinto de cero en esta región. El efecto es sutil, pero es un sesgo de todos modos. (La fidelidad deRbfes incluso mejor con la opción predeterminada de funciones radiales, denominadas multiquadratic).

Función maligna y muestreo superior

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

Están comenzando a aparecer claras diferencias entre los tres métodos. Al interp2dobservar las gráficas de superficie, hay claros extremos falsos que aparecen en la salida de (observe las dos jorobas en el lado derecho de la superficie trazada). Si bien griddatay Rbfparecen producir resultados similares a primera vista, este último parece producir un mínimo más profundo cerca [0.4, -0.4]que está ausente de la función subyacente.

Sin embargo, hay un aspecto crucial en el que Rbfes 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 griddatarompe la simetría de los puntos de muestra, que ya es débilmente visible en el caso suave.

Funcionamiento fluido y datos dispersos

La mayoría de las veces se desea realizar una interpolación en datos dispersos. Por eso espero que estas pruebas sean más importantes. Como se muestra arriba, los puntos de la 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, para empezar, debe considerar si tiene sentido interpolar sus datos sin procesar.

Salida para la función suave:

Ahora ya hay un pequeño espectáculo de terror. Recorté la salida de interp2da entre [-1, 1]exclusivamente para trazar, 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 enormes regiones ruidosas donde el método se descompone por completo. El segundo caso de griddatareproduce la forma bastante bien, pero observe las regiones blancas en el borde del gráfico de contorno. Esto se debe al hecho de que griddatasolo 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. ² Teniendo en cuenta estas características, Rbfparece funcionar mejor.

Función maligna y datos dispersos

Y el momento que todos hemos estado esperando:

No es una gran sorpresa que se interp2drinda. De hecho, durante la llamada interp2ddeberías esperar que algunos amigos se RuntimeWarningquejen de la imposibilidad de construir el spline. En cuanto a los otros dos métodos, Rbfparece producir el mejor resultado, incluso cerca de los límites del dominio donde se extrapola el resultado.

Así que permítanme decir algunas palabras sobre los tres métodos, en orden decreciente de preferencia (para que el peor sea el que tenga menos probabilidades de ser leído por alguien).

scipy.interpolate.Rbf

La Rbfclase significa "funciones de base radial". Para ser honesto, nunca consideré 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 splines (ver más adelante), el uso viene en dos pasos: primero uno crea una Rbfinstancia de clase 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 ascendente 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_rbftiene la misma forma x_densey y_densesin ningún esfuerzo. También tenga en cuenta que Rbfadmite dimensiones arbitrarias para la interpolación.

Entonces, scipy.interpolate.Rbf

• produce una salida con buen comportamiento incluso para datos de entrada locos
• admite la 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 en absoluto)
• crea un interpolador como primer paso, por lo que evaluarlo en varios puntos de salida supone menos esfuerzo adicional
• puede tener puntos de salida de forma arbitraria (a diferencia de estar restringido a mallas rectangulares, ver más adelante)
• propenso a preservar la simetría de los datos de entrada
• soporta múltiples tipos de funciones radiales para la palabra clave function: multiquadric, inverse, gaussian, linear, cubic, quintic, thin_platey arbitraria definida por el usuario

scipy.interpolate.griddata

Mi antiguo favorito, griddataes un caballo de batalla general para la interpolación en dimensiones arbitrarias. No realiza una 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 torpe. Los puntos de entrada deben especificarse en una matriz de formas [N, D]en Ddimensiones. 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 zdatos de entrada también deben aplanarse. 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 helpde griddataes engañoso, 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 una palabra, scipy.interpolate.griddata

• produce una salida con buen comportamiento incluso para datos de entrada locos
• admite la interpolación en dimensiones superiores
• no realiza extrapolación, se puede establecer un valor único 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
• admite interpolación lineal y de vecino más cercano en dimensiones arbitrarias, cúbica en 1d y 2d. Uso de interpolación lineal y de vecino más cercano NearestNDInterpolatory LinearNDInterpolatorbajo el capó, respectivamente. La interpolación cúbica 1d utiliza un spline, la interpolación cúbica 2d se utiliza CloughTocher2DInterpolatorpara 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 hablando de interp2dsus parientes 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 porque se usa específicamente para la interpolación bidimensional, pero sospecho que este es, con mucho, el caso más común de interpolación multivariante.

En lo que respecta a la sintaxis, interp2des similar Rbfen que primero necesita construir una instancia de interpolación, que se puede 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 fueran 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 usarlo interp2des poner todas las mallas 2d en la llamada de interpolación, lo que lleva a un consumo de memoria explosivo y, con suerte, a un apresurado MemoryError.

Ahora, el mayor problema con interp2des que a menudo no funciona. Para entender esto, tenemos que mirar debajo del capó. Resulta que interp2des 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 0.17.0 de scipy) hay un buen comentariointerpolate/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 bisplrephay 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 subyacentes interp2dno están destinadas realmente a realizar interpolaciones. Pueden ser suficientes para datos con un comportamiento suficientemente bueno, pero en circunstancias realistas probablemente querrá usar algo más.

Solo para concluir, interpolate.interp2d

• puede provocar artefactos incluso con datos bien temperados
• es específicamente para problemas bivariados (aunque existe el límite interpnpara 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 supone 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
• admite interpolación lineal, cúbica y quíntica
• podría violar la simetría de los datos de entrada

¹ Estoy bastante seguro de que las funciones básicas cubicy el lineartipo de Rbfno se corresponden exactamente con los otros interpoladores del mismo nombre.
² Estos NaN también son la razón por la que el gráfico de superficie parece tan extraño: históricamente, matplotlib ha tenido dificultades para trazar objetos 3D complejos con la 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 al frente. Este es un problema de visualización y no de interpolación.