python - tutorial - que es numpy

Interpolación lineal Python 4D en una cuadrícula rectangular

Question

Necesito interpolar linealmente los datos de temperatura en 4 dimensiones (latitud, longitud, altitud y tiempo).
El número de puntos es bastante alto (360x720x50x8) y necesito un método rápido para calcular la temperatura en cualquier punto del espacio y el tiempo dentro de los límites de datos.

He intentado usar scipy.interpolate.LinearNDInterpolator pero utilizar Qhull para la triangulación es ineficiente en una cuadrícula rectangular y lleva horas completarlo.

Al leer este boleto SciPy , la solución parecía estar implementando un nuevo interpolador nd usando la interp1d estándar para calcular una mayor cantidad de puntos de datos, y luego usar un enfoque de "vecino más cercano" con el nuevo conjunto de datos.

Esto, sin embargo, lleva mucho tiempo nuevamente (minutos).

¿Hay una forma rápida de interpolar datos en una cuadrícula rectangular en 4 dimensiones sin que tome minutos lograrlo?

Pensé en usar interp1d 4 veces sin calcular una mayor densidad de puntos, pero dejándolo para que el usuario llame con las coordenadas, pero no entiendo cómo hacerlo.

De lo contrario, ¿sería una opción escribir mi propio interpolador 4D específico para mis necesidades?

Aquí está el código que he estado usando para probar esto:

Usando scipy.interpolate.LinearNDInterpolator :

import numpy as np
from scipy.interpolate import LinearNDInterpolator
lats = np.arange(-90,90.5,0.5)
lons = np.arange(-180,180,0.5)
alts = np.arange(1,1000,21.717)
time = np.arange(8)
data = np.random.rand(len(lats)*len(lons)*len(alts)*len(time)).reshape((len(lats),len(lons),len(alts),len(time)))
coords = np.zeros((len(lats),len(lons),len(alts),len(time),4))
coords[...,0] = lats.reshape((len(lats),1,1,1))
coords[...,1] = lons.reshape((1,len(lons),1,1))
coords[...,2] = alts.reshape((1,1,len(alts),1))
coords[...,3] = time.reshape((1,1,1,len(time)))
coords = coords.reshape((data.size,4))
interpolatedData = LinearNDInterpolator(coords,data)

Usando scipy.interpolate.interp1d :

import numpy as np
from scipy.interpolate import LinearNDInterpolator
lats = np.arange(-90,90.5,0.5)
lons = np.arange(-180,180,0.5)
alts = np.arange(1,1000,21.717)
time = np.arange(8)
data = np.random.rand(len(lats)*len(lons)*len(alts)*len(time)).reshape((len(lats),len(lons),len(alts),len(time)))
interpolatedData = np.array([None, None, None, None])
interpolatedData[0] = interp1d(lats,data,axis=0)
interpolatedData[1] = interp1d(lons,data,axis=1)
interpolatedData[2] = interp1d(alts,data,axis=2)
interpolatedData[3] = interp1d(time,data,axis=3)

¡Muchas gracias por su ayuda!

Answer 1

En el mismo ticket que ha vinculado, hay una implementación de ejemplo de lo que ellos llaman interpolación de producto tensor , que muestra la forma correcta de anidar llamadas recursivas a interp1d . Esto es equivalente a la interpolación quadrilinear si elige el parámetro kind=''linear'' predeterminado para sus interp1d ''s.

Si bien esto puede ser lo suficientemente bueno, esto no es una interpolación lineal, y habrá términos de orden superior en la función de interpolación, ya que esta imagen de la entrada de la wikipedia en la interpolación bilineal muestra:

Esto bien puede ser lo suficientemente bueno para lo que está buscando, pero hay aplicaciones en las que se prefiere una interpolación lineal triangulada, realmente a trozos. Si realmente lo necesita, existe una manera fácil de solucionar la lentitud de qhull.

Una vez que se ha configurado LinearNDInterpolator , hay dos pasos para obtener un valor interpolado para un punto dado:

1. averiguar dentro de qué triángulo (hipertetraedro 4D en su caso) el punto es, y
2. interpolar usando las coordenadas baricéntricas del punto relativo a los vértices como pesos.

Probablemente no quiera meterse con las coordenadas baricéntricas, así que mejor déjelo en LinearNDInterpolator . Pero sabes algunas cosas sobre la triangulación. Sobre todo eso, debido a que tienes una cuadrícula regular, dentro de cada hipercubo la triangulación va a ser la misma. Entonces, para interpolar un solo valor, primero se podría determinar en qué subcubo se encuentra su punto, construir un LinearNDInterpolator con los 16 vértices de ese cubo, y usarlo para interpolar su valor:

from itertools import product
def interpolator(coords, data, point) :
dims = len(point)
indices = []
sub_coords = []
for j in xrange(dims) :
idx = np.digitize([point[j]], coords[j])[0]
indices += [[idx - 1, idx]]
sub_coords += [coords[j][indices[-1]]]
indices = np.array([j for j in product(*indices)])
sub_coords = np.array([j for j in product(*sub_coords)])
sub_data = data[list(np.swapaxes(indices, 0, 1))]
li = LinearNDInterpolator(sub_coords, sub_data)
return li([point])[0]
>>> point = np.array([12.3,-4.2, 500.5, 2.5])
>>> interpolator((lats, lons, alts, time), data, point)
0.386082399091

Esto no puede funcionar en datos vectorizados, ya que eso requeriría almacenar un LinearNDInterpolator para cada posible subcubo, y aunque probablemente sería más rápido que triangular todo, aún sería muy lento.

Answer 2

Para cosas muy similares, utilizo Scientific.Functions.Interpolation.InterpolatingFunction .

 import numpy as np
from Scientific.Functions.Interpolation import InterpolatingFunction
lats = np.arange(-90,90.5,0.5)
lons = np.arange(-180,180,0.5)
alts = np.arange(1,1000,21.717)
time = np.arange(8)
data = np.random.rand(len(lats)*len(lons)*len(alts)*len(time)).reshape((len(lats),len(lons),len(alts),len(time)))
axes = (lats, lons, alts, time)
f = InterpolatingFunction(axes, data)

Ahora puede dejar que el usuario llame a InterpolatingFunction con coordenadas:

>>> f(0,0,10,3)
0.7085675631375401

InterpolatingFunction tiene buenas características adicionales, como integración y división.

Sin embargo, no estoy seguro de si la interpolación es lineal. Tendría que buscar en la fuente del módulo para averiguarlo.

Answer 3

scipy.ndimage.map_coordinates es un buen interpolador rápido para grillas uniformes (todas las casillas del mismo tamaño). Ver multivariate-spline-interpolation-in-python-scipy en SO para una descripción clara.

Para cuadrículas rectangulares no uniformes, una envoltura simple Intergrid mapea / escala de forma no uniforme a grillas uniformes, luego hace map_coordinates. En un caso de prueba 4d como el suyo, se necesitan aproximadamente 1 μseg por consulta:

Intergrid: 1000000 points in a (361, 720, 47, 8) grid took 652 msec