Interpolación sobre una cuadrícula irregular

Question

Entonces, tengo tres matrices numpy que almacenan latitud, longitud y algún valor de propiedad en una cuadrícula, es decir, tengo LAT (y, x), LON (y, x) y , diga la temperatura T (y, x), para algunos límites de x y y. La cuadrícula no es necesariamente regular, de hecho, es tripolar.

Quiero interpolar estos valores de propiedad (temperatura) en un grupo de diferentes puntos lat/lon (almacenados como lat1 (t), lon1 (t), por aproximadamente 10,000 t ...) que no caen en los puntos de cuadrícula reales. 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 forma hábil, preferiblemente rápida de hacer esto? No puedo evitar pensar que la respuesta es obvia ... ¡Gracias!

Answer 1

Tratar la combinación de ponderación inversa de la distancia y scipy.spatial.KDTree descrito en SO inverse-distance-weighted-idw-interpolation-with-python. Kd-trees funciona bien en 2d 3D ..., la ponderación de distancia inversa es uniforme y local, y el número k = de vecinos más cercanos se puede variar a la velocidad/precisión de compromiso.

Answer 2

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.

Answer 3

¿Estoy en lo cierto al pensar que sus cuadrículas de datos se ven algo como esto (el rojo es el dato antiguo, el azul es el nuevo dato interpolado)?

alt text http://www.geekops.co.uk/photos/0000-00-02%20%28Forum%20images%29/DataSeparation.png

Esto podría ser un enfoque ligeramente por fuerza bruta-ish, pero ¿qué pasa con la prestación de su información existente como un mapa de bits (OpenGL va a hacer una simple interpolación de colores para usted con las opciones correctas configurados y que podrían hacer 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.

Answer 4

Hay un montón de opciones aquí, cuál es el mejor dependerá de sus datos ... Sin embargo no sé de una solución fuera de la caja para usted

Usted dice que su entrada los datos provienen de datos tripolares. Hay tres casos principales sobre cómo se pueden estructurar estos datos.

1. Muestreado desde una cuadrícula 3D en el espacio tripolar, proyectado de nuevo a datos 2d LAT, LON.
2. Muestreado desde una cuadrícula 2d en espacio tripolar, proyectado en 2d datos LAT LON.
3. Los datos no estructurados en el espacio tripolar proyectadas en datos LON LAT 2d

La forma más fácil de estos es 2. En lugar de la interpolación en el espacio LON LAT, "sólo" transformar su punto nuevo en el espacio de origen e interpolar 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 volver a hacer 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, utilizaría una estrategia no estructurada, ya que me preocuparía tener que gestionar la búsqueda entre celdas con proyecciones superpuestas. Elegir la celda "correcta" sería difícil.

Answer 5

Hay un nice inverse distance example by Roger Veciana i Rovira junto con un código que usa GDAL para escribir a geotiff si te interesa.

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í, al mismo tiempo tenga cuidado con qué proyección se utiliza para sus datos.

Una copia de su algoritmo con la prueba de:

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() 

Answer 6

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

partir de la descripción:

BIVAR es una biblioteca FORTRAN90 el cual interpola datos de dos variables dispersas, 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 ser evaluado en otros puntos en el plano.