Inverse Distance Weighted (IDW) interpolazione con Python

Question

la domanda: Qual è il modo migliore per calcolare la distanza inversa ponderate (IDW) interpolazione in Python, per le posizioni dei punti?

Alcuni background: Attualmente sto usando RPy2 per interfacciarlo con R e il suo modulo gstat. Sfortunatamente, il modulo gstat è in conflitto con arcgisscripting che ho ottenuto eseguendo l'analisi basata su RPy2 in un processo separato. Anche se questo problema è stato risolto in una versione recente/futura, e l'efficienza può essere migliorata, vorrei comunque rimuovere la mia dipendenza dall'installazione di R.

Il sito Web gstat fornisce un eseguibile autonomo, che è più facile da pacchetto con il mio script python, ma mi auguro ancora una soluzione Python che non richieda più scritture su disco e lanci processi esterni. Il numero di chiamate alla funzione di interpolazione, di insiemi separati di punti e valori, può avvicinarsi a 20.000 nell'elaborazione che sto eseguendo.

Ho specificamente bisogno di interpolare per punti, quindi usare la funzione IDW in ArcGIS per generare suoni di raster è ancora peggio che usare R, in termini di prestazioni ..... a meno che non ci sia un modo per mascherare in modo efficiente solo i punti Ho bisogno. Anche con questa modifica, non mi aspetterei che le prestazioni siano così eccezionali. Vedrò questa opzione come un'altra alternativa. AGGIORNARE: Il problema qui è che sei legato alla dimensione della cella che stai utilizzando. Se riduci la dimensione della cella per ottenere una maggiore precisione, l'elaborazione richiede molto tempo. È inoltre necessario eseguire il follow-up estraendo per punti ..... un metodo troppo brutto se si desidera valori per punti specifici.

Ho visto lo scipy documentation, ma non sembra che ci sia un modo semplice per calcolare IDW.

Sto pensando di implementare la mia implementazione, possibilmente utilizzando alcune delle funzionalità di scipy per individuare i punti più vicini e calcolare le distanze.

Mi manca qualcosa di ovvio? C'è un modulo Python che non ho visto che fa esattamente quello che voglio? Creare la mia implementazione con l'aiuto di scipy è una scelta saggia?

Answer 1

cambiati 20 ottobre: questa classe Invdisttree combina la ponderazione della distanza inversa e scipy.spatial.KDTree.
Dimentica la risposta originale a forza bruta; questo è imho il metodo di scelta per l'interpolazione di dati sparsi.

""" invdisttree.py: inverse-distance-weighted interpolation using KDTree 
    fast, solid, local 
""" 
from __future__ import division 
import numpy as np 
from scipy.spatial import cKDTree as KDTree 
    # http://docs.scipy.org/doc/scipy/reference/spatial.html 

__date__ = "2010-11-09 Nov" # weights, doc 

#............................................................................... 
class Invdisttree: 
    """ inverse-distance-weighted interpolation using KDTree: 
invdisttree = Invdisttree(X, z) -- data points, values 
interpol = invdisttree(q, nnear=3, eps=0, p=1, weights=None, stat=0) 
    interpolates z from the 3 points nearest each query point q; 
    For example, interpol[ a query point q ] 
    finds the 3 data points nearest q, at distances d1 d2 d3 
    and returns the IDW average of the values z1 z2 z3 
     (z1/d1 + z2/d2 + z3/d3) 
     /(1/d1 + 1/d2 + 1/d3) 
     = .55 z1 + .27 z2 + .18 z3 for distances 1 2 3 

    q may be one point, or a batch of points. 
    eps: approximate nearest, dist <= (1 + eps) * true nearest 
    p: use 1/distance**p 
    weights: optional multipliers for 1/distance**p, of the same shape as q 
    stat: accumulate wsum, wn for average weights 

How many nearest neighbors should one take ? 
a) start with 8 11 14 .. 28 in 2d 3d 4d .. 10d; see Wendel's formula 
b) make 3 runs with nnear= e.g. 6 8 10, and look at the results -- 
    |interpol 6 - interpol 8| etc., or |f - interpol*| if you have f(q). 
    I find that runtimes don't increase much at all with nnear -- ymmv. 

p=1, p=2 ? 
    p=2 weights nearer points more, farther points less. 
    In 2d, the circles around query points have areas ~ distance**2, 
    so p=2 is inverse-area weighting. For example, 
     (z1/area1 + z2/area2 + z3/area3) 
     /(1/area1 + 1/area2 + 1/area3) 
     = .74 z1 + .18 z2 + .08 z3 for distances 1 2 3 
    Similarly, in 3d, p=3 is inverse-volume weighting. 

Scaling: 
    if different X coordinates measure different things, Euclidean distance 
    can be way off. For example, if X0 is in the range 0 to 1 
    but X1 0 to 1000, the X1 distances will swamp X0; 
    rescale the data, i.e. make X0.std() ~= X1.std() . 

A nice property of IDW is that it's scale-free around query points: 
if I have values z1 z2 z3 from 3 points at distances d1 d2 d3, 
the IDW average 
    (z1/d1 + z2/d2 + z3/d3) 
    /(1/d1 + 1/d2 + 1/d3) 
is the same for distances 1 2 3, or 10 20 30 -- only the ratios matter. 
In contrast, the commonly-used Gaussian kernel exp(- (distance/h)**2) 
is exceedingly sensitive to distance and to h. 

    """ 
# anykernel(dj/av dj) is also scale-free 
# error analysis, |f(x) - idw(x)| ? todo: regular grid, nnear ndim+1, 2*ndim 

    def __init__(self, X, z, leafsize=10, stat=0): 
     assert len(X) == len(z), "len(X) %d != len(z) %d" % (len(X), len(z)) 
     self.tree = KDTree(X, leafsize=leafsize) # build the tree 
     self.z = z 
     self.stat = stat 
     self.wn = 0 
     self.wsum = None; 

    def __call__(self, q, nnear=6, eps=0, p=1, weights=None): 
      # nnear nearest neighbours of each query point -- 
     q = np.asarray(q) 
     qdim = q.ndim 
     if qdim == 1: 
      q = np.array([q]) 
     if self.wsum is None: 
      self.wsum = np.zeros(nnear) 

     self.distances, self.ix = self.tree.query(q, k=nnear, eps=eps) 
     interpol = np.zeros((len(self.distances),) + np.shape(self.z[0])) 
     jinterpol = 0 
     for dist, ix in zip(self.distances, self.ix): 
      if nnear == 1: 
       wz = self.z[ix] 
      elif dist[0] < 1e-10: 
       wz = self.z[ix[0]] 
      else: # weight z s by 1/dist -- 
       w = 1/dist**p 
       if weights is not None: 
        w *= weights[ix] # >= 0 
       w /= np.sum(w) 
       wz = np.dot(w, self.z[ix]) 
       if self.stat: 
        self.wn += 1 
        self.wsum += w 
      interpol[jinterpol] = wz 
      jinterpol += 1 
     return interpol if qdim > 1 else interpol[0] 

#............................................................................... 
if __name__ == "__main__": 
    import sys 

    N = 10000 
    Ndim = 2 
    Nask = N # N Nask 1e5: 24 sec 2d, 27 sec 3d on mac g4 ppc 
    Nnear = 8 # 8 2d, 11 3d => 5 % chance one-sided -- Wendel, mathoverflow.com 
    leafsize = 10 
    eps = .1 # approximate nearest, dist <= (1 + eps) * true nearest 
    p = 1 # weights ~ 1/distance**p 
    cycle = .25 
    seed = 1 

    exec "\n".join(sys.argv[1:]) # python this.py N= ... 
    np.random.seed(seed) 
    np.set_printoptions(3, threshold=100, suppress=True) # .3f 

    print "\nInvdisttree: N %d Ndim %d Nask %d Nnear %d leafsize %d eps %.2g p %.2g" % (
     N, Ndim, Nask, Nnear, leafsize, eps, p) 

    def terrain(x): 
     """ ~ rolling hills """ 
     return np.sin((2*np.pi/cycle) * np.mean(x, axis=-1)) 

    known = np.random.uniform(size=(N,Ndim)) ** .5 # 1/(p+1): density x^p 
    z = terrain(known) 
    ask = np.random.uniform(size=(Nask,Ndim)) 

#............................................................................... 
    invdisttree = Invdisttree(known, z, leafsize=leafsize, stat=1) 
    interpol = invdisttree(ask, nnear=Nnear, eps=eps, p=p) 

    print "average distances to nearest points: %s" % \ 
     np.mean(invdisttree.distances, axis=0) 
    print "average weights: %s" % (invdisttree.wsum/invdisttree.wn) 
     # see Wikipedia Zipf's law 
    err = np.abs(terrain(ask) - interpol) 
    print "average |terrain() - interpolated|: %.2g" % np.mean(err) 

    # print "interpolate a single point: %.2g" % \ 
    #  invdisttree(known[0], nnear=Nnear, eps=eps) 

Answer 2

Edit: @Denis è giusto, un Rbf lineare (ad esempio scipy.interpolate.Rbf con "la funzione = 'lineare'") non è la stessa di IDW ...

(Nota, tutti questi useranno quantità eccessive di memoria, se si utilizza un gran numero di punti)

Ecco un semplice exampe di IDW:

def simple_idw(x, y, z, xi, yi): 
    dist = distance_matrix(x,y, xi,yi) 

    # In IDW, weights are 1/distance 
    weights = 1.0/dist 

    # Make weights sum to one 
    weights /= weights.sum(axis=0) 

    # Multiply the weights for each interpolated point by all observed Z-values 
    zi = np.dot(weights.T, z) 
    return zi 

considerando che, ecco quello che una Rbf lineare sarebbe:

def linear_rbf(x, y, z, xi, yi): 
    dist = distance_matrix(x,y, xi,yi) 

    # Mutual pariwise distances between observations 
    internal_dist = distance_matrix(x,y, x,y) 

    # Now solve for the weights such that mistfit at the observations is minimized 
    weights = np.linalg.solve(internal_dist, z) 

    # Multiply the weights for each interpolated point by the distances 
    zi = np.dot(dist.T, weights) 
    return zi 

(Usando la funzione distance_matrix qui :)

def distance_matrix(x0, y0, x1, y1): 
    obs = np.vstack((x0, y0)).T 
    interp = np.vstack((x1, y1)).T 

    # Make a distance matrix between pairwise observations 
    # Note: from <http://stackoverflow.com/questions/1871536> 
    # (Yay for ufuncs!) 
    d0 = np.subtract.outer(obs[:,0], interp[:,0]) 
    d1 = np.subtract.outer(obs[:,1], interp[:,1]) 

    return np.hypot(d0, d1) 

Mettere tutto insieme in un bel esempio di copia-incolla cede alcuni appezzamenti Confronto veloce: Homemade IDW example plot http://www.geology.wisc.edu/~jkington/homemade_idw.pngHomemade linear RBF example plot http://www.geology.wisc.edu/~jkington/homemade_rbf.pngScipy's linear RBF example plot http://www.geology.wisc.edu/~jkington/scipy_rbf.png

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.interpolate import Rbf 

def main(): 
    # Setup: Generate data... 
    n = 10 
    nx, ny = 50, 50 
    x, y, z = map(np.random.random, [n, n, n]) 
    xi = np.linspace(x.min(), x.max(), nx) 
    yi = np.linspace(y.min(), y.max(), ny) 
    xi, yi = np.meshgrid(xi, yi) 
    xi, yi = xi.flatten(), yi.flatten() 

    # Calculate IDW 
    grid1 = simple_idw(x,y,z,xi,yi) 
    grid1 = grid1.reshape((ny, nx)) 

    # Calculate scipy's RBF 
    grid2 = scipy_idw(x,y,z,xi,yi) 
    grid2 = grid2.reshape((ny, nx)) 

    grid3 = linear_rbf(x,y,z,xi,yi) 
    print grid3.shape 
    grid3 = grid3.reshape((ny, nx)) 


    # Comparisons... 
    plot(x,y,z,grid1) 
    plt.title('Homemade IDW') 

    plot(x,y,z,grid2) 
    plt.title("Scipy's Rbf with function=linear") 

    plot(x,y,z,grid3) 
    plt.title('Homemade linear Rbf') 

    plt.show() 

def simple_idw(x, y, z, xi, yi): 
    dist = distance_matrix(x,y, xi,yi) 

    # In IDW, weights are 1/distance 
    weights = 1.0/dist 

    # Make weights sum to one 
    weights /= weights.sum(axis=0) 

    # Multiply the weights for each interpolated point by all observed Z-values 
    zi = np.dot(weights.T, z) 
    return zi 

def linear_rbf(x, y, z, xi, yi): 
    dist = distance_matrix(x,y, xi,yi) 

    # Mutual pariwise distances between observations 
    internal_dist = distance_matrix(x,y, x,y) 

    # Now solve for the weights such that mistfit at the observations is minimized 
    weights = np.linalg.solve(internal_dist, z) 

    # Multiply the weights for each interpolated point by the distances 
    zi = np.dot(dist.T, weights) 
    return zi 


def scipy_idw(x, y, z, xi, yi): 
    interp = Rbf(x, y, z, function='linear') 
    return interp(xi, yi) 

def distance_matrix(x0, y0, x1, y1): 
    obs = np.vstack((x0, y0)).T 
    interp = np.vstack((x1, y1)).T 

    # Make a distance matrix between pairwise observations 
    # Note: from <http://stackoverflow.com/questions/1871536> 
    # (Yay for ufuncs!) 
    d0 = np.subtract.outer(obs[:,0], interp[:,0]) 
    d1 = np.subtract.outer(obs[:,1], interp[:,1]) 

    return np.hypot(d0, d1) 


def plot(x,y,z,grid): 
    plt.figure() 
    plt.imshow(grid, extent=(x.min(), x.max(), y.max(), y.min())) 
    plt.hold(True) 
    plt.scatter(x,y,c=z) 
    plt.colorbar() 

if __name__ == '__main__': 
    main()