[GIS] What’s wrong with the IDW interpolation function (Python)

algorithminterpolationinverse-distance-weightedpython

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.

Best Answer

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
Related Question