Inverse Distance Weighted (IDW) Interpolation with Python
changed 20 Oct: this class Invdisttree combines inverse-distance weighting and
scipy.spatial.KDTree.
Forget the original brute-force answer;
this is imho the method of choice for scattered-data interpolation.
""" 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 )
Python geospatial interpolation (meteorological data)
Found solution here:
Inverse Distance Weighted (IDW) Interpolation with Python
IDW interpolation is more than enough in my case, but @user6386471, thanks for your contribution!
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
(Using the distance_matrix function here:)
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)
Inverse Distance Weighted interpolation in a grid(Python)
I figured out how to solve it.
First i calculated the distance between each labeled cell and each unlabeled cell, using this function
def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * asin(sqrt(a))
# Radius of earth in kilometers is 6371
km = 6371* c
return km
Then i created a function called iwd that apply the inverse weighted distance interpolation using the standard formula.
After that i was able to obtain a a label for each cell that result in the following image:
How to do R spatial interpolation with an inverse distance weighting error
Even though this error was caused by a typo, (formual = instead of formula = in the call to idw), the error message is a little opaque, so in case anyone comes across a similar error in the future, I thought I would show what it means and why it comes about.
In R you can set up a new generic that allows you to call different methods according to the class of the object that is passed to a function (plot is the classic example here, where you can call on the same function to plot different types of object, but completely different code is called depending on the object type in order to generate a sensible plot).
We can set up a boring example of a generic. Suppose we want to have a generic function called boring that, when passed two numbers, will simply add them, but when passed two character strings, will paste them together into a single string.
setGeneric("boring", function(x, y, ...) standardGeneric("boring"))
#> [1] "boring"
setMethod("boring", c("numeric", "numeric"), function(x, y) x + y)
setMethod("boring", c("character", "character"), function(x, y) paste(x, y))
This works as expected:
boring(x = 1, y = 2)
#> [1] 3
boring(x = "a", y = "b")
#> [1] "a b"
Because we have specified the type of arguments that are allowed, if we try to pass in one number and one character, we get an error telling us that there is no method available for that signature:
boring(x = 1, y = "b")
#> Error in (function (classes, fdef, mtable) : unable to find an inherited method
#> for function 'boring' for signature '"numeric", "character"'
However, if we completely miss out x when calling boring (say by accidentally using z instead of x due to a typo), we don't get the standard R error telling us that x is missing with no default value. Instead, we get the error message reported in this question:
boring(z = 1, y = 2)
#> Error in (function (classes, fdef, mtable) : unable to find an inherited method
#> for function 'boring' for signature '"missing", "numeric"'
So the error just means you have missed out or misnamed a parameter to a generic function.
Created on 2021-11-04 by the reprex package (v2.0.0)
Creat a inverse distance weighting Function in VBA
As far as I can see, all you want to do is set values to be zero if they are empty:
Function IDWW(Value1, Value2, Value3, Dist1,Dist2,Dist3)
Dim a1 As Variant
Dim b1 As Variant
Dim a2 As Variant
Dim b2 As Variant
Dim a3 As Variant
Dim b3 As Variant
If Value1 <> "" And Dist1 <> "" Then
a1 = Value1 / (Dist1) ^ 2
b1 = 1 / (Dist1) ^ 2
Else
a1 = 0
b1 = 0
End If
If Value2 <> "" And Dist2 <> "" Then
a2 = Value2 / (Dist2) ^ 2
b2 = 1 / (Dist2) ^ 2
Else
a2 = 0
b2 = 0
End If
If Value3 <> "" And Dist3 <> "" Then
a3 = Value3 / (Dist3) ^ 2
b3 = 1 / (Dist3) ^ 2
Else
a3 = 0
b3 = 0
End If
'Avoid a problem if all 3 distances are empty
If b1 + b2 + b3 = 0 Then
IDWW = 0
Else
IDWW = (a1+a2+a3) / (b1+b2+b3)
End If
End Function
Related Topics
Multithreaded Blas in Python/Numpy
Read from a Log File as It's Being Written Using Python
Check If String Has Date, Any Format
Python - Rolling Functions for Groupby Object
Python: Call a Function from String Name
How to Add Items into a Numpy Array
What's a Faster Operation, Re.Match/Search or Str.Find
The Problem with Installing Pil Using Virtualenv or Buildout
Meaning of Inter_Op_Parallelism_Threads and Intra_Op_Parallelism_Threads
Concatenate a List of Pandas Dataframes Together
Use and Meaning of "In" in an If Statement
Clicking on Svg Using Selenium Python
Creating Over 20 Unique Legend Colors Using Matplotlib
Moving X-Axis to the Top of a Plot in Matplotlib
How to Print a List in Python "Nicely"
Run Child Processes as Different User from a Long Running Python Process