What's wrong with my IDW interpolation function (Python)?

Original Question

Can anyone see what's wrong with my custom IDW interpolation below?

I know there are existing libraries that can do this for me, but for various reasons including educational and stubbornness, I'm trying to create a pure Python Inverse Distance Weighted Interpolation function. I'm basing myself off the "basic form" of the Shepard algorithm from this Wiki article:

and

The way I interpret these equations is `weightofknownpoint = 1/(knowntounknowndistance to the power of the P parameter)` and `IDWforunknownpoint = sumforallknownpoints(weightofknownpoint*valueofknownpoint/thesumof(allknownweights))`

The Python code I came up with produces interpolated values that are way off (the values spike in between known points instead of doing a smooth transition). To limit the length of this post I'm not going to show all of my code, so you won't be able to run this on your own. Just assume that everything else is correct (btw I know the code might not be very optimized but I'll deal with that later); is there anything about the code below that doesn't look right in terms of what the IDW algorithm is supposed to be doing, that might explain my odd results:

``````#ASSUMPTIONS
#math module is imported
#unknowncell and knowncell refer to custom created instances with a .x and .y position properties, and a .value property for their data value to be interpolated
#the unknowncell argument that I use for the _CalcIDWvalue function are all the unknown cell instances in the grid
#knowncells = list of known cell instances (ie the known points used to inform the interpolated value)
#sensitivity = 3 #ie the power parameter

def CalcIDWvalue(unknowncell):
calclist = []
for knowncell in knowncells:
if cell.x == knowncell.x and cell.y == knowncell.y:
continue
#prep
def _getweight(unknowncell,knowncell):
dist = math.sqrt((unknowncell.x-knowncell.x)**2 + (unknowncell.y-knowncell.y)**2)
return 1/float(dist**sensitivity)
weight = _getweight(unknowncell,knowncell)
sumofallweights = sum([_getweight(unknowncell,tempknowncell) for tempknowncell in knowncells
if unknowncell.x != tempknowncell.x and unknowncell.y != tempknowncell.y])
#calculate
tempcalc = (weight*knowncell.value) / sumofallweights
calclist.append(tempcalc)
return sum(calclist)
``````

My input grid is this:

``````[0, None, 1, None, 2, None, 3, None, 4, None, 5]
[None, None, None, None, None, None, None, None, None, None, None]
[0, None, 1, None, 2, None, 3, None, 4, None, 5]
[None, None, None, None, None, None, None, None, None, None, None]
[0, None, 1, None, 2, None, 3, None, 4, None, 5]
[None, None, None, None, None, None, None, None, None, None, None]
[0, None, 1, None, 2, None, 3, None, 4, None, 5]
[None, None, None, None, None, None, None, None, None, None, None]
[0, None, 1, None, 2, None, 3, None, 4, None, 5]
[None, None, None, None, None, None, None, None, None, None, None]
[0, None, 1, None, 2, None, 3, None, 4, None, 5]
``````

And my incorrectly interpolated output grid is:

``````[0.0, 4.92, 1.0, 11.76, 2.0, 18.79, 3.0, 26.47, 4.0, 36.33, 5.0]
[1.42, 0.66, 5.54, 1.56, 9.92, 2.5, 14.64, 3.44, 20.43, 4.34, 39.83]
[0.0, 3.17, 1.0, 7.51, 2.0, 12.0, 3.0, 16.82, 4.0, 22.8, 5.0]
[1.43, 0.67, 5.22, 1.57, 9.28, 2.5, 13.68, 3.43, 19.11, 4.33, 36.8]
[0.0, 3.03, 1.0, 7.1, 2.0, 11.32, 3.0, 15.86, 4.0, 21.53, 5.0]
[1.44, 0.68, 5.17, 1.57, 9.15, 2.5, 13.48, 3.43, 18.84, 4.32, 36.14]
[0.0, 3.03, 1.0, 7.1, 2.0, 11.32, 3.0, 15.86, 4.0, 21.53, 5.0]
[1.43, 0.67, 5.22, 1.57, 9.28, 2.5, 13.68, 3.43, 19.11, 4.33, 36.8]
[0.0, 3.17, 1.0, 7.51, 2.0, 12.0, 3.0, 16.82, 4.0, 22.8, 5.0]
[1.42, 0.66, 5.54, 1.56, 9.92, 2.5, 14.64, 3.44, 20.43, 4.34, 39.83]
[0.0, 4.92, 1.0, 11.76, 2.0, 18.79, 3.0, 26.47, 4.0, 36.33, 5.0]
``````

Update

I found out what was wrong and posted the solution as an answer below. Another answer by Rickhg23s has further helped optimize the speed of the code, so see his answer for the final most up-to-date version of the IDW code.

• For your input grid, should we assume that the distance between each grid point in a row/column is 1? Commented Mar 22, 2014 at 19:10
• Yes, but that's only because that was the particular test example I was working with. The regularly spaced grid is really just regularly spaced points "placed onto" an empty grid. The IDW function should be able to deal with any types of locations of the points/knowncells. Anyway, I found the solution so check out my answer below :) Commented Mar 23, 2014 at 5:04

In case you want to make it go faster by reducing duplicated computations:

``````def CalcIDWvalue(unknowncell, knowncells):
weighted_values_sum = 0.0
sum_of_weights = 0.0
neg_half_sens = -sensitivity/2.0
for knowncell in knowncells:
weight = ((unknowncell.x-knowncell.x)**2 + (unknowncell.y-knowncell.y)**2)**neg_half_sens
sum_of_weights += weight
weighted_values_sum += weight * knowncell.value
return weighted_values_sum / sum_of_weights
``````
• Nicely spotted, appreciate the extra pair of eyes. It did speed up with about 10%, not much but better than nothing, and besides the function was only used once so was a bit redundant. I took the liberty to change your version of the code to what you suggested, pending your approval, since that is the accepted answer. I had previously come across some of the IDW posts you linked to, but I had only seen ones relying on Numpy and/or Scipy. I'll see if I can dig some more later maybe, but for now I'm happy with the function. For huge datasets, one could simply limit knowncells by N nearest/dist. Commented Mar 24, 2014 at 23:40
• Hey, why quit now? I streamlined it a bit more. Commented Mar 25, 2014 at 4:09
• That's nuts. Your most recent streamline didn't seem that radical, but it appears it gave another speed increase of about 30%. In practical terms, all your're optimizations have now decreased the time for interpolating 20 points on a 600*600 grid down from 40 secs to 13 secs. That makes quite a bit of difference, with better support for interpolating larger sized grids and sample points. You're a wiz at this game! Commented Mar 25, 2014 at 13:42
• Python is supposed to be pretty efficient. In other (less efficient) systems, interpolating 20 points via IDW onto a 600 by 600 grid will be almost instantaneous. 13 seconds is still ridiculously slow. I don't know Python well enough to suggest what is wrong, but the first place I would look for a speedup is the algorithm, which ought to loop over the data points at the highest level (at least when there are so few of them), and next in terms of employing vectorized operations. With this approach, `R` does the interpolation in 0.3 second on my workstation (30 seconds for a 6000 by 6000 grid). Commented Mar 25, 2014 at 18:13
• @whuber I've found R's vector operations to be very fast indeed. And you're right, to get maximum performance, minimizing/vectorizing calculations with good algorithms is critical. Commented Mar 25, 2014 at 22:04

Awesome, well I found out which parts of the code were wrong and got it to work. I ended up removing two segments in the code that only did some action when the unknowncell and the knowncell were different from each other, which is pointless to test for since that should never be the case. Although obviously it must have had some effect since it managed to mess with the IDW values. I might not understand why, but at least I ended up fixing the bug.

Thus with this pure-Python function for calculating the IDW value, and some other helper functions and code for dealing with grids and putting the point data onto a grid that I'm not showing here, I was able to create the following point-to-grid interpolation (the image itself offcourse required the use of Numpy and PIL):

Note that for the interpolation seen in the image I used a much larger grid input than in the original question, and the points were randomly scattered instead of evenly spaced. The speed was reasonable, taking about 40 seconds for up to 600*600 grids (improved to 20 seconds with the accepted answer by rickhg12hs), but for 1000*1000 it started to struggle and take much much longer, so it's certainly not the most effective interpolation alternative, but for small tasks it will do. Here is the final working IDW function that did it:

Update: Note that the code below is not the preferred one. See instead the accepted answer by Rickhg12hs.

``````#ASSUMPTIONS
#math module is imported
#unknowncell and knowncell refer to custom created instances with a .x and .y position properties, and a .value property for their data value to be interpolated
#the unknowncell argument that I use for the _CalcIDWvalue function are all the unknown cell instances in the grid
#knowncells = list of known cell instances (ie the known points used to inform the interpolated value)
#sensitivity = 3 #ie the power parameter

def CalcIDWvalue(unknowncell):
calclist = []
for knowncell in knowncells:
#prep
def _getweight(unknowncell,knowncell):
dist = math.sqrt((unknowncell.x-knowncell.x)**2 + (unknowncell.y-knowncell.y)**2)
return 1/float(dist**sensitivity)
weight = _getweight(unknowncell,knowncell)
sumofallweights = sum([_getweight(unknowncell,tempknowncell) for tempknowncell in knowncells])
#calculate
tempcalc = (weight*knowncell.value) / sumofallweights
calclist.append(tempcalc)
return sum(calclist)
``````