python - guia - Interpolación rápida de datos de cuadrícula
qgis manual (2)
Tengo una gran np.ndarray 3d de datos que representa una variable física muestreada sobre un volumen en una forma de cuadrícula regular (como en el valor de la matriz [0,0,0] representa el valor en coordenadas físicas (0,0,0) )).
Me gustaría ir a un espaciado de cuadrícula más fino mediante la interpolación de los datos en la grilla aproximada. En este momento estoy usando interpolación lineal de malla cuadrada, pero es bastante lenta (~ 90 segundos para una matriz de 20x20x20). Es un poco sobredimensionado para mis propósitos, permitiendo el muestreo aleatorio de los datos de volumen. ¿Hay algo que pueda aprovechar mis datos regularmente espaciados y el hecho de que solo hay un conjunto limitado de puntos específicos a los que quiero interpolar?
Gran respuesta de Joe. De acuerdo con su sugerencia, creé el paquete regulargrid ( https://pypi.python.org/pypi/regulargrid/ , fuente en https://github.com/JohannesBuchner/regulargrid )
Proporciona soporte para grillas cartesianas n-dimensionales (según sea necesario aquí) a través de la muy rápida scipy.ndimage.map_coordinates para escalas de coordenadas arbitrarias.
¡Por supuesto! Hay dos opciones que hacen cosas diferentes, pero ambas explotan la naturaleza de cuadrícula regular de los datos originales.
El primero es scipy.ndimage.zoom . Si solo quieres producir una grilla regular más densa basada en la interpolación de los datos originales, este es el camino a seguir.
El segundo es scipy.ndimage.map_coordinates . Si desea interpolar algunos (o muchos) puntos arbitrarios en sus datos, pero aún así explotar la naturaleza de cuadrícula regular de los datos originales (por ejemplo, no se requiere quadtree), es el camino a seguir.
"Acercar" una matriz ( scipy.ndimage.zoom )
Como un ejemplo rápido (Esto usará interpolación cúbica. Use order=1 para bilinear, order=0 para lo más cercano, etc.):
import numpy as np
import scipy.ndimage as ndimage
data = np.arange(9).reshape(3,3)
print ''Original:/n'', data
print ''Zoomed by 2x:/n'', ndimage.zoom(data, 2)
Esto produce:
Original:
[[0 1 2]
[3 4 5]
[6 7 8]]
Zoomed by 2x:
[[0 0 1 1 2 2]
[1 1 1 2 2 3]
[2 2 3 3 4 4]
[4 4 5 5 6 6]
[5 6 6 7 7 7]
[6 6 7 7 8 8]]
Esto también funciona para matrices 3D (y nD). Sin embargo, tenga en cuenta que si hace zoom por 2x, por ejemplo, hará zoom en todos los ejes.
data = np.arange(27).reshape(3,3,3)
print ''Original:/n'', data
print ''Zoomed by 2x gives an array of shape:'', ndimage.zoom(data, 2).shape
Esto produce:
Original:
[[[ 0 1 2]
[ 3 4 5]
[ 6 7 8]]
[[ 9 10 11]
[12 13 14]
[15 16 17]]
[[18 19 20]
[21 22 23]
[24 25 26]]]
Zoomed by 2x gives an array of shape: (6, 6, 6)
Si tiene algo así como una imagen RGB de 3 bandas que desea ampliar, puede hacerlo especificando una secuencia de tuplas como factor de zoom:
print ''Zoomed by 2x along the last two axes:''
print ndimage.zoom(data, (1, 2, 2))
Esto produce:
Zoomed by 2x along the last two axes:
[[[ 0 0 1 1 2 2]
[ 1 1 1 2 2 3]
[ 2 2 3 3 4 4]
[ 4 4 5 5 6 6]
[ 5 6 6 7 7 7]
[ 6 6 7 7 8 8]]
[[ 9 9 10 10 11 11]
[10 10 10 11 11 12]
[11 11 12 12 13 13]
[13 13 14 14 15 15]
[14 15 15 16 16 16]
[15 15 16 16 17 17]]
[[18 18 19 19 20 20]
[19 19 19 20 20 21]
[20 20 21 21 22 22]
[22 22 23 23 24 24]
[23 24 24 25 25 25]
[24 24 25 25 26 26]]]
Interpolación arbitraria de datos de cuadrícula regular utilizando map_coordinates
Lo primero que hay que entender sobre map_coordinates es que funciona en coordenadas de píxel (por ejemplo, igual que map_coordinates la matriz, pero los valores pueden ser flotantes). Según su descripción, esto es exactamente lo que quiere, pero si a menudo confunde a las personas. Por ejemplo, si tiene coordenadas x, y, z del "mundo real", deberá transformarlas en coordenadas de "píxeles" basadas en índices.
En cualquier caso, digamos que queríamos interpolar el valor en la matriz original en la posición 1.2, 0.3, 1.4.
Si está pensando en esto en términos del caso de imagen RGB anterior, la primera coordenada corresponde a la "banda", la segunda a la "fila" y la última a la "columna". El orden que corresponde a lo que depende completamente de cómo decida estructurar sus datos, pero voy a usarlos como coordenadas "z, y, x", ya que hace que la comparación con la matriz impresa sea más fácil de visualizar.
import numpy as np
import scipy.ndimage as ndimage
data = np.arange(27).reshape(3,3,3)
print ''Original:/n'', data
print ''Sampled at 1.2, 0.3, 1.4:''
print ndimage.map_coordinates(data, [[1.2], [0.3], [1.4]])
Esto produce:
Original:
[[[ 0 1 2]
[ 3 4 5]
[ 6 7 8]]
[[ 9 10 11]
[12 13 14]
[15 16 17]]
[[18 19 20]
[21 22 23]
[24 25 26]]]
Sampled at 1.2, 0.3, 1.4:
[14]
Una vez más, esto es interpolación cúbica por defecto. Utilice la order kwarg para controlar el tipo de interpolación.
Vale la pena señalar aquí que todas las operaciones de scipy.ndimage conservan el tipo de la matriz original. Si desea obtener resultados en coma flotante, deberá convertir la matriz original como flotante:
In [74]: ndimage.map_coordinates(data.astype(float), [[1.2], [0.3], [1.4]])
Out[74]: array([ 13.5965])
Otra cosa que puede notar es que el formato de coordenadas interpoladas es bastante engorroso para un solo punto (por ejemplo, espera una matriz 3xN en lugar de una matriz Nx3). Sin embargo, es posiblemente más agradable cuando tienes secuencias de coordenadas. Por ejemplo, considere el caso del muestreo a lo largo de una línea que pasa por el "cubo" de datos:
xi = np.linspace(0, 2, 10)
yi = 0.8 * xi
zi = 1.2 * xi
print ndimage.map_coordinates(data, [zi, yi, xi])
Esto produce:
[ 0 1 4 8 12 17 21 24 0 0]
Este también es un buen lugar para mencionar cómo se manejan las condiciones de contorno. Por defecto, cualquier cosa que esté fuera de la matriz se establece en 0. Por lo tanto, los dos últimos valores de la secuencia son 0 . (es decir, zi es> 2 para los dos últimos elementos).
Si queremos que los puntos fuera de la matriz sean, digamos -999 (No podemos usar nan ya que esta es una matriz de enteros. Si quieres nan , tendrás que convertir a flotantes):
In [75]: ndimage.map_coordinates(data, [zi, yi, xi], cval=-999)
Out[75]: array([ 0, 1, 4, 8, 12, 17, 21, 24, -999, -999])
Si quisiéramos que devolviera el valor más cercano para los puntos fuera del conjunto, lo haríamos:
In [76]: ndimage.map_coordinates(data, [zi, yi, xi], mode=''nearest'')
Out[76]: array([ 0, 1, 4, 8, 12, 17, 21, 24, 25, 25])
También puede usar "reflect" y "wrap" como modos de contorno, además de "nearest" y el "constant" predeterminado. Estos son bastante autoexplicativos, pero intente experimentar un poco si está confundido.
Por ejemplo, vamos a interpolar una línea a lo largo de la primera fila de la primera banda en la matriz que se extiende por el doble de la distancia de la matriz:
xi = np.linspace(0, 5, 10)
yi, zi = np.zeros_like(xi), np.zeros_like(xi)
El valor predeterminado es:
In [77]: ndimage.map_coordinates(data, [zi, yi, xi])
Out[77]: array([0, 0, 1, 2, 0, 0, 0, 0, 0, 0])
Compare esto con:
In [78]: ndimage.map_coordinates(data, [zi, yi, xi], mode=''reflect'')
Out[78]: array([0, 0, 1, 2, 2, 1, 2, 1, 0, 0])
In [78]: ndimage.map_coordinates(data, [zi, yi, xi], mode=''wrap'')
Out[78]: array([0, 0, 1, 2, 0, 1, 1, 2, 0, 1])
¡Espero que eso aclare un poco las cosas!