Interpolación sobre una rejilla irregular.

Entonces, tengo tres matrices numpy que almacenan la latitud, la longitud y algún valor de propiedad en una cuadrícula, es decir, tengo LAT (y, x), LON (y, x) y, digamos, 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 propiedades (temperatura) en un grupo de diferentes puntos lat / lon (almacenados como lat1 (t), lon1 (t), para aproximadamente 10,000 t …) que no caen en los puntos de la cuadrícula reales . He intentado matplotlib.mlab.griddata, pero eso lleva demasiado tiempo (no está realmente diseñado para lo que estoy haciendo, después de todo). También probé scipy.interpolate.interp2d, pero obtengo un MemoryError (mis grids son aproximadamente 400×400).

¿Hay algún tipo de manera hábil, preferiblemente rápida de hacer esto? No puedo evitar pensar que la respuesta es algo obvio … ¡Gracias!

Pruebe la combinación de ponderación de distancia inversa y scipy.spatial.KDTree descrita en SO inversa-distancia-ponderada-idw-interpolación-con-python . Los árboles-kd funcionan bien en 2d 3d …, la ponderación en distancias inversas es suave y local, y el k = número de vecinos más cercanos se puede variar a la velocidad / precisión del intercambio.

Hay un buen ejemplo de distancia inversa de Roger Veciana i Rovira junto con algún código que usa GDAL para escribir al geotiff si te gusta eso.

Esto es de una cuadrícula regular, pero suponiendo que proyecte los datos primero en una cuadrícula de píxeles con pyproj o algo así, tenga cuidado con la proyección que se usa 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 un montón de opciones aquí, la mejor dependerá de sus datos … Sin embargo, no conozco una solución lista para usar

Usted dice que sus datos de entrada son de datos tripolares. Hay tres casos principales de cómo se podrían estructurar estos datos.

  1. Muestreo de una cuadrícula 3d en el espacio tripolar, proyectado de nuevo a datos LAT, LON 2d.
  2. Muestreo de una cuadrícula 2d en el espacio tripolar, proyectado en datos LON LON 2d.
  3. Datos no estructurados en el espacio tripolar proyectados en datos LAT LON 2d

El más fácil de estos es 2. En lugar de interpolar en el espacio LAT LON, “simplemente” transforme su punto de nuevo en el espacio de origen e interpole allí.

Otra opción que funciona para 1 y 2 es buscar las celdas que se asignan desde el espacio tripolar para cubrir su punto de muestra. (Puede utilizar una estructura de tipo BSP o de cuadrícula para acelerar esta búsqueda) Elija una de las celdas e interpole dentro de ella.

Finalmente, hay un montón de opciones de interpolación no estructuradas … pero tienden a ser lentas. Uno de mis favoritos personales es usar una interpolación lineal de los N puntos más cercanos, encontrando que esos N puntos se pueden hacer nuevamente con cuadrícula o un BSP. Otra buena opción es alinear los puntos desestructurados de Delauney e interpolarlos 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 celdas con proyecciones superpuestas. Elegir la celda “correcta” sería difícil.

Le sugiero que eche un vistazo a las características de interpolación de GRASS (un paquete de SIG de código abierto) ( http://grass.ibiblio.org/gdp/html_grass62/v.surf.bspline.html ). No está en Python, pero puede volver a implementarlo o interactuar con el código C.

¿Estoy en lo cierto al pensar que las cuadrículas de datos tienen un aspecto similar a este (el rojo es el antiguo, el azul es el nuevo interpolado)?

texto alternativo http://sofes.miximages.com/python/DataSeparation.png

Esto podría ser un enfoque ligeramente brusco, pero ¿qué hay de representar sus datos existentes como un bitmap? (Opengl hará una simple interpolación de colores con las opciones correctas configuradas y puede representar los datos como triangularjs, que deberían ser bastante rápidos ). A continuación, puede muestrear píxeles en las ubicaciones de los nuevos puntos.

Alternativamente, puede ordenar su primer conjunto de puntos espacialmente y luego encontrar los puntos antiguos más cercanos que rodean su nuevo punto e interpolar según las distancias a esos puntos.

Hay una biblioteca de FORTRAN llamada BIVAR , que es muy adecuada para este problema. Con algunas modificaciones, puede hacer que sea utilizable en python usando f2py.

De la descripción:

BIVAR es una biblioteca de FORTRAN90 que interpola datos bivariados dispersos, de Hiroshi Akima.

BIVAR acepta un conjunto de (X, Y) puntos de datos dispersos en 2D, con valores de datos Z asociados, y puede construir una función de interpolación suave Z (X, Y), que concuerda con los datos dados, y puede ser evaluada en Otros puntos en el plano.