python - transferencia - librerias de scipy
Interpolación sobre una cuadrÃcula irregular (6)
¿Estoy en lo cierto al pensar que sus cuadrículas de datos se parecen a esto (el rojo es el dato antiguo, el azul es el nuevo dato interpolado)?
Este podría ser un enfoque un poco brute-force-ish, pero ¿qué pasa con renderizar sus datos existentes como un mapa de bits? (OpenGL hará una interpolación simple de colores para usted con las opciones correctas configuradas y podría representar los datos como triángulos que deberían ser bastante rápidos ) A continuación, puede muestrear píxeles en las ubicaciones de los nuevos puntos.
Alternativamente, podría ordenar su primer conjunto de puntos espacialmente y luego encontrar los puntos antiguos más cercanos a su nuevo punto e interpolar en función de las distancias a esos puntos.
Entonces, tengo tres matrices numpy que almacenan latitud, longitud y algún valor de propiedad en una grilla, es decir, tengo LAT (y, x), LON (y, x) y, por ejemplo, temperatura T (y, x ), para algunos límites de x y y. La cuadrícula no es necesariamente regular, de hecho, es tripolar.
Luego quiero interpolar estos valores de propiedad (temperatura) en un grupo de diferentes puntos lat / lon (almacenados como lat1 (t), lon1 (t), por alrededor de 10,000 t ...) que no caen en los puntos reales de la grilla . Intenté con matplotlib.mlab.griddata, pero eso lleva demasiado tiempo (en realidad, no está diseñado para lo que estoy haciendo, después de todo). También probé scipy.interpolate.interp2d, pero obtengo un MemoryError (mis cuadrículas son aproximadamente 400x400).
¿Hay algún tipo de hábil, preferiblemente rápida forma de hacer esto? No puedo evitar pensar que la respuesta es obvia ... ¡Gracias!
Aquí hay un montón de opciones, cuál es la mejor dependerá de sus datos ... Sin embargo, no sé de una solución lista para usar.
Usted dice que sus datos de entrada provienen de datos tripolares. Hay tres casos principales sobre cómo se pueden estructurar estos datos.
- Tomado de una grilla 3d en el espacio tripolar, proyectado de nuevo a 2d LAT, datos LON.
- Tomado de una grilla 2d en espacio tripolar, proyectado en 2d LAT LON data.
- Datos no estructurados en el espacio tripolar proyectados en datos 2d LAT LON
El más fácil de estos es 2. En lugar de interpolar en el espacio LAT LON, "solo" transforme su punto nuevamente en el espacio de origen e interpole allí.
Otra opción que funciona para 1 y 2 es buscar las celdas que se correlacionan desde el espacio tripolar para cubrir su punto de muestra. (Puede usar una estructura de tipo BSP o de cuadrícula para acelerar esta búsqueda) Elija una de las celdas e interpola dentro de ella.
Finalmente, hay un montón de opciones de interpolación no estructuradas ... pero tienden a ser lentas. Uno de mis favoritos es utilizar una interpolación lineal de los N puntos más cercanos, encontrando que los N puntos se pueden hacer nuevamente con grillado o un BSP. Otra buena opción es triangular Delauney los puntos no estructurados e interpolar en la malla triangular resultante.
Personalmente, si mi malla fuera el caso 1, usaría una estrategia no estructurada ya que me preocuparía tener que manejar la búsqueda a través de las células con proyecciones superpuestas. Elegir la celda "correcta" sería difícil.
Hay un buen ejemplo de distancia inversa de Roger Veciana i Rovira junto con un código que usa GDAL para escribir a geotiff si te interesa eso.
Esto es de gruesa a una cuadrícula regular, pero suponiendo que proyecte los datos primero a una cuadrícula de píxeles con pyproj o algo así, todo el tiempo teniendo cuidado con qué proyección se utiliza para sus datos.
Una copia de su algoritmo con prueba :
from math import pow
from math import sqrt
import numpy as np
import matplotlib.pyplot as plt
def pointValue(x,y,power,smoothing,xv,yv,values):
nominator=0
denominator=0
for i in range(0,len(values)):
dist = sqrt((x-xv[i])*(x-xv[i])+(y-yv[i])*(y-yv[i])+smoothing*smoothing);
#If the point is really close to one of the data points, return the data point value to avoid singularities
if(dist<0.0000000001):
return values[i]
nominator=nominator+(values[i]/pow(dist,power))
denominator=denominator+(1/pow(dist,power))
#Return NODATA if the denominator is zero
if denominator > 0:
value = nominator/denominator
else:
value = -9999
return value
def invDist(xv,yv,values,xsize=100,ysize=100,power=2,smoothing=0):
valuesGrid = np.zeros((ysize,xsize))
for x in range(0,xsize):
for y in range(0,ysize):
valuesGrid[y][x] = pointValue(x,y,power,smoothing,xv,yv,values)
return valuesGrid
if __name__ == "__main__":
power=1
smoothing=20
#Creating some data, with each coodinate and the values stored in separated lists
xv = [10,60,40,70,10,50,20,70,30,60]
yv = [10,20,30,30,40,50,60,70,80,90]
values = [1,2,2,3,4,6,7,7,8,10]
#Creating the output grid (100x100, in the example)
ti = np.linspace(0, 100, 100)
XI, YI = np.meshgrid(ti, ti)
#Creating the interpolation function and populating the output matrix value
ZI = invDist(xv,yv,values,100,100,power,smoothing)
# Plotting the result
n = plt.normalize(0.0, 100.0)
plt.subplot(1, 1, 1)
plt.pcolor(XI, YI, ZI)
plt.scatter(xv, yv, 100, values)
plt.title(''Inv dist interpolation - power: '' + str(power) + '' smoothing: '' + str(smoothing))
plt.xlim(0, 100)
plt.ylim(0, 100)
plt.colorbar()
plt.show()
Hay una biblioteca de FORTRAN llamada BIVAR , que es muy adecuada para este problema. Con algunas modificaciones, puede hacerlo utilizable en python usando f2py.
De la descripción:
BIVAR es una biblioteca FORTRAN90 que interpola datos bivariados dispersos, por Hiroshi Akima.
BIVAR acepta un conjunto de puntos de datos (X, Y) dispersos en 2D, con valores de datos Z asociados, y es capaz de construir una función de interpolación suave Z (X, Y), que concuerda con los datos dados, y puede evaluarse a otros puntos en el plano.
Le sugiero que eche un vistazo a las funciones de interpolación GRASS (un paquete de GIS de código abierto) ( http://grass.ibiblio.org/gdp/html_grass62/v.surf.bspline.html ). No está en Python, pero puedes volver a implementarlo o interactuar con el código C.
Pruebe la combinación de ponderación de distancia inversa y scipy.spatial.KDTree descrita en SO inverse-distance-weighted-idw-interpolation-with-python . Kd-trees KD funcionan muy bien en 2d 3D ..., la ponderación de distancia inversa es suave y local, y el número de k = de vecinos más cercanos puede variarse a la velocidad / precisión de compensación.