Interpolación con distancia ponderada inversa (IDW) con Python

La pregunta: ¿Cuál es la mejor manera de calcular la interpolación de distancia ponderada (IDW) en Python, para ubicaciones de puntos?

Algunos antecedentes: Actualmente estoy usando RPy2 para interactuar con R y su módulo gstat. Desafortunadamente, el módulo gstat está en conflicto con arcgisscripting, que resolví ejecutando el análisis basado en RPy2 en un proceso separado. Incluso si este problema se resuelve en una versión reciente / futura, y la eficiencia se puede mejorar, todavía me gustaría eliminar mi dependencia de la instalación de R.

El sitio web gstat proporciona un ejecutable independiente, que es más fácil de empaquetar con mi script de Python, pero todavía espero una solución de Python que no requiera múltiples escrituras en el disco y que inicie procesos externos. El número de llamadas a la función de interpolación, de conjuntos separados de puntos y valores, puede aproximarse a 20,000 en el procesamiento que estoy realizando.

Necesito específicamente interpolar puntos, por lo que usar la función IDW en ArcGIS para generar rásteres suena aún peor que usar R, en términos de rendimiento … a menos que haya una manera de enmascarar de manera eficiente solo los puntos que necesito. Incluso con esta modificación, no esperaría que el rendimiento fuera tan bueno. Veré esta opción como otra alternativa. ACTUALIZACIÓN: El problema aquí es que está vinculado al tamaño de celda que está utilizando. Si reduce el tamaño de la celda para obtener una mayor precisión, el procesamiento toma mucho tiempo. También debe hacer un seguimiento extrayendo puntos … sobre todo un método feo si desea valores para puntos específicos.

He mirado la documentación de scipy , pero no parece que haya una manera directa de calcular IDW.

Estoy pensando en rodar mi propia implementación, posiblemente utilizando alguna de las funciones de scipy para ubicar los puntos más cercanos y calcular distancias.

¿Me estoy perdiendo algo obvio? ¿Hay un módulo de python que no haya visto que haga exactamente lo que quiero? ¿Crear una implementación propia con la ayuda de scipy es una elección inteligente?

cambiado el 20 de octubre: esta clase Invdisttree combina la ponderación de distancia inversa y scipy.spatial.KDTree .
Olvida la respuesta original de la fuerza bruta; Este es el método de elección para la interpolación de datos dispersos.

""" 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= eg 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, ie 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 zs 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 ) 

Edición: @Denis tiene razón, un Rbf lineal (por ejemplo, scipy.interpolate.Rbf con “function = ‘linear'”) no es lo mismo que IDW …

(Tenga en cuenta que todos estos usarán cantidades excesivas de memoria si está usando una gran cantidad de puntos).

Aquí hay un ejemplo simple de 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 que, esto es lo que sería un Rbf lineal:

 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 función distance_matrix aquí 🙂

 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  # (Yay for ufuncs!) d0 = np.subtract.outer(obs[:,0], interp[:,0]) d1 = np.subtract.outer(obs[:,1], interp[:,1]) return np.hypot(d0, d1) 

Poniéndolo todo junto en un bonito ejemplo de copiar y pegar, se obtienen algunas plots de comparación rápida: Parcela de ejemplo IDW hecha en casa http://www.geology.wisc.edu/~jkington/homemade_idw.png Parcela de ejemplo RBF lineal hecha en casa http: // www. geology.wisc.edu/~jkington/homemade_rbf.png gráfico de ejemplo de RBF lineal de Scipy 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  # (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()