Interpolacja nad nieregularną siatką

Istnieje ładny przykład odwróconej odległości Roger Veciana i Rovira wraz z kodem używającym GDAL do pisania do geotiffa, jeśli to lubisz.

Jest to gruba do zwykłej siatki, ale zakładając, że najpierw rzutujesz dane do siatki pikseli za pomocą pyproj lub czegoś takiego, cały czas uważając, jaka projekcja jest używana dla Twoich danych.

Kopia jego algorytmu z testem :

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

Jest tu kilka opcji, która z nich jest najlepsza, zależy od Twoich danych... Jednak nie znam gotowego rozwiązania dla Ciebie

Mówisz, że dane wejściowe pochodzą z danych tripolarnych. Istnieją trzy główne przypadki dotyczące struktury tych danych.

1. próbkowane z siatki 3d w przestrzeni trójpolarnej, rzutowane z powrotem do danych 2D LAT, LON.
2. próbkowane z siatki 2d w przestrzeni trójpolarnej, rzutowane na dane 2D LAT LON.
3. Dane nieustrukturyzowane w przestrzeni trójpolarnej projektowane na dane 2d lat lon

Najprostszym z nich jest 2. Zamiast interpolować w przestrzeni LAT LON," po prostu " przekształć swój punkt z powrotem w Przestrzeń źródłową i tam Interpoluj.

Inną opcją, która działa dla 1 i 2, jest wyszukiwanie komórek, które mapują z przestrzeni trójpolarnej, aby pokryć punkt próbki. (Możesz użyć struktury typu BSP lub siatki, aby przyspieszyć to wyszukiwanie) wybierz jedną z komórek i Interpoluj wewnątrz niej.

Wreszcie jest sterta niestrukturalnych opcje interpolacji .. ale są powolne. Moim ulubionym jest użycie liniowej interpolacji najbliższych N punktów, znalezienie tych N punktów można ponownie zrobić za pomocą gridding lub BSP. Inną dobrą opcją jest Delauney triangulacji punktów niestrukturalnych i interpolacji na powstałej siatki trójkątne.

Osobiście, gdyby moja siatka była przypadkiem 1, użyłbym niestrukturalnej strategii, ponieważ martwiłbym się o przeszukiwanie komórek z nakładającymi się projekcjami. Wybór "właściwej" komórki byłby trudny.

Proponuję rzucić okiem na Grass (open source GIS package) interpolation features ( http://grass.ibiblio.org/gdp/html_grass62/v.surf.bspline.html ). nie jest w Pythonie, ale można go ponownie zaimplementować lub interfejs z kodem C.

Czy mam rację myśląc, że Twoje siatki danych wyglądają tak (czerwony to stare dane, niebieski to nowe dane interpolowane)?

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

Może to być nieco brutalne podejście, ale co z renderowaniem istniejących danych jako bitmapy (opengl zrobi prostą interpolację kolorów dla Ciebie z odpowiednimi opcjami skonfigurowanymi i możesz renderować dane jako Trójkąty który powinien być dość szybki). Następnie można próbkować piksele w miejscach nowych punktów.

Alternatywnie można posortować pierwszy zbiór punktów przestrzennie, a następnie znaleźć najbliższe stare punkty otaczające nowy punkt i interpolować na podstawie odległości do tych punktów.

Istnieje Biblioteka FORTRAN o nazwieBIVAR , która jest bardzo odpowiednia dla tego problemu. Dzięki kilku modyfikacjom możesz uczynić go użytecznym w Pythonie za pomocą f2py.

Z opisu:

BIVAR jest biblioteką FORTRAN90, która interpoluje rozproszone dane biwarialne, autorstwa Hiroshiego Akimy.

BIVAR przyjmuje zbiór (X, Y) punktów danych rozproszonych w 2D, z powiązanymi wartościami danych Z i jest w stanie skonstruować płynną funkcję interpolacyjną Z (X, Y), która zgadza się z podanymi danymi i mogą być oceniane w innych punktach na płaszczyźnie.