Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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:

enter image description here

and

enter image description here

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.

share|improve this question
    
For your input grid, should we assume that the distance between each grid point in a row/column is 1? –  rickhg12hs Mar 22 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 :) –  Karim Bahgat Mar 23 at 5:04

2 Answers 2

up vote 3 down vote accepted

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
share|improve this answer
1  
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. –  Karim Bahgat Mar 24 at 23:40
1  
Hey, why quit now? I streamlined it a bit more. –  rickhg12hs Mar 25 at 4:09
1  
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! –  Karim Bahgat Mar 25 at 13:42
1  
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). –  whuber Mar 25 at 18:13
1  
@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. –  rickhg12hs Mar 25 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):

enter image description here

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)
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.