python polynomiale Interpolation sur une grille irrégulière



interpolation python (5)

Donc, j'ai trois tableaux numpy qui stockent la latitude, la longitude et une valeur de propriété sur une grille - c'est-à-dire que j'ai LAT (y, x), LON (y, x) et, disons, température T (y, x ), pour certaines limites de x et y. La grille n'est pas forcément régulière - en fait, elle est tripolaire.

Je souhaite ensuite interpoler ces valeurs de propriété (température) sur un groupe de points de lat / lon différents (stockés sous la forme lat1 (t), lon1 (t), pour environ 10 000 t ...) qui ne correspondent pas aux points de grille réels . J'ai essayé matplotlib.mlab.griddata, mais cela prend beaucoup trop de temps (ce n'est pas vraiment conçu pour ce que je fais, après tout). J'ai aussi essayé scipy.interpolate.interp2d, mais je reçois une MemoryError (mes grilles font environ 400x400).

Y a-t-il une sorte de slick, de préférence un moyen rapide de le faire? Je ne peux pas m'empêcher de penser que la réponse est évidente… Merci !!



Answer #2

Ai-je raison de penser que vos grilles de données ressemblent à ceci (les anciennes données sont en rouge, les nouvelles données interpolées en bleu)?

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

Cela peut être une approche légèrement brutale, mais qu'en est-il du rendu de vos données existantes en tant que bitmap (opengl effectuera une interpolation simple des couleurs avec les bonnes options configurées et vous pourrez rendre les données sous forme de triangles assez rapides) ). Vous pouvez ensuite échantillonner des pixels aux emplacements des nouveaux points.

Vous pouvez également trier votre premier ensemble de points spatialement, puis rechercher les anciens points les plus proches de votre nouveau point et les interpoler en fonction des distances qui les séparent.


Answer #3

Roger Veciana i Rovira propose un exemple de distance inverse avec du code utilisant GDAL pour écrire à geotiff si vous y êtes.

Ceci est grossier par rapport à une grille régulière, mais en supposant que vous projetez d'abord les données sur une grille de pixels avec pyproj ou autre, tout en faisant attention à la projection utilisée pour vos données.

Une copie de son algorithme avec test :

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 #4

Il existe une bibliothèque FORTRAN appelée BIVAR , qui convient parfaitement à ce problème. Avec quelques modifications, vous pouvez le rendre utilisable en python en utilisant f2py.

De la description:

BIVAR est une bibliothèque FORTRAN90 qui interpole des données bivariées dispersées, par Hiroshi Akima.

BIVAR accepte un ensemble de (X, Y) points de données dispersés en 2D, avec des valeurs de données Z associées, et est capable de construire une fonction d'interpolation lisse Z (X, Y), qui correspond aux données données et peut être évaluée à d'autres points dans l'avion.


Answer #5

Il y a un tas d'options ici, le meilleur dépendra de vos données ... Cependant, je ne connais pas de solution prête à l'emploi pour vous

Vous dites que vos données d'entrée proviennent de données tripolaires. Il existe trois principaux cas de structuration de ces données.

  1. Echantillonné à partir d'une grille 3D dans l'espace tripolaire, projeté vers les données 2D LAT, LON.
  2. Echantillonné à partir d'une grille 2d dans l'espace tripolaire, projeté en données 2D LAT LON.
  3. Données non structurées dans l'espace tripolaire projetées dans des données LAT LON 2D

Le plus simple d'entre eux est 2. Au lieu d'interpoler dans l'espace LAT LON, "juste" transformez votre point dans l'espace source et y interpolez.

Une autre option qui fonctionne pour 1 et 2 consiste à rechercher les cellules correspondant à l’espace tripolaire pour couvrir votre point d’échantillon. (Vous pouvez utiliser une structure de type BSP ou grille pour accélérer cette recherche) Sélectionnez l'une des cellules et interpolez-la.

Enfin, il y a un tas d'options d'interpolation non structurées .. mais elles ont tendance à être lentes. Un de mes favoris personnels est d'utiliser une interpolation linéaire des N points les plus proches, trouver ces N points peut encore être fait avec un maillage ou un BSP. Une autre bonne option consiste à Delauney trianguler les points non structurés et à interpoler sur le maillage triangulaire résultant.

Personnellement, si mon maillage était le cas 1, j'utiliserais une stratégie non structurée car je m'inquiéterais de devoir gérer la recherche à travers des cellules avec des projections qui se chevauchent. Choisir la "bonne" cellule serait difficile.





interpolation