I have been looking for python script to interpolate my data. I found the code that can do my question.

I have a couple of questions in modifying the code below. 1) My data is ordanized an a way as given by data.txt. It has lat, long and values. I want the code read my data and process from my text data. 2) I also have the lat and long grid points. I want to do the interpolation from ths grid point data and append the interpolated values with the grid point data provided. 3) I have NoData/NULL values in mt given data, hence I want the code to do the interpolation only from the stations that have data only.

data.xtx

``XV YV v1 v2 v2 v4 v5 v6 v710 10 1 2 4 3 NA 410 15 4 4 3 NA NA NA30 35 NA NA 1 NA 5 185 20 4 NA 4 3 10 NA15 15 NA 5 4 NA NA 525 10 7 8 7 5 10 NA``

My grid point data looks

grid.txt

`` X Y10 1010 1510 2010 2510 3010 3510 4015 1015 1515 2015 2515 3015 3515 4020 1020 1520 2020 2520 3020 3520 4025 1025 1525 2025 2525 3025 3525 4030 1030 1530 2030 2530 3030 3530 4035 1035 1535 2035 2535 3035 3535 4040 1040 1540 2040 2540 3040 3540 40``

The code I am trying is

``#! /usr/bin/pythonfrom math import powfrom math import sqrtimport numpy as npdef pointValue(x,y,power,smoothing,xv,yv,values):nominator=0denominator=0for i in range(0,len(values)):dist = sqrt((x-xv[i])*(x-xv[i])+(y-yv[i])*(y-yv[i])+smoothing*smoothing);#If the point is really close to one of the data points, return the data point value to avoid singularitiesif(dist<0.0000000001):return values[i]nominator=nominator+(values[i]/pow(dist,power))denominator=denominator+(1/pow(dist,power))#Return NODATA if the denominator is zeroif denominator > 0:value = nominator/denominatorelse:value = -9999return valuedef invDist(xv,yv,values,xsize=100,ysize=100,power=2,smoothing=0):valuesGrid = np.zeros((ysize,xsize))for x in range(0,xsize):for y in range(0,ysize):valuesGrid[y][x] = pointValue(x,y,power,smoothing,xv,yv,values)return valuesGridif __name__ == "__main__":power=1smoothing=20#Creating some data, with each coodinate and the values stored in separated listsxv = [10,60,40,70,10,50,20,70,30,60]yv = [10,20,30,30,40,50,60,70,80,90]values = [1,2,2,3,4,6,7,7,8,10]#Creating the output grid (100x100, in the example)ti = np.linspace(0, 100, 100)XI, YI = np.meshgrid(ti, ti)#Creating the interpolation function and populating the output matrix valueZI = invDist(xv,yv,values,100,100,power,smoothing)print ZI``

Are you positively sure you want IDW? IDW has the nasty property that when you interpolate something, you will have to take all points into account. It is also quite bad at interpolating things in many cases. So, if you just need a decent interpolation for your data, please have a look at `scipy.interpolate.griddata` It might be exactly what you need.

On the other hand, if you want to interpolate a list of unknown points by using IDW, it can be vectorised to be rather short (but remember that this is still O(m*n), where m and n are the numbers of unknown and known points).

``````import numpy as np

# inverse power to use
pow = 2

# some points:
known_pts = np.array([[0,0], [1,0], [3,0], [4,0],
[0,1], [2,1], [3,1], [4,1],
[1,3], [2,3], [3,3], [4,3]])
# and some data associated with them
known_values = known_pts[:,0] + known_pts[:,1]

# unknown points
unknown_pts = np.array([[2,0], [1,1],  [0,2], [1,2], [2,2], [3,2], [4,2], [0,3]])

# calculate all distances from unknown points to known points:
# (the array will have as many rows as there are unknown points and as many columns
#  as there are known points)
dist = np.sqrt((known_pts[:,0][None,:]-unknown_pts[:,0][:,None])**2 + (known_pts[:,1][None,:]-unknown_pts[:,1][:,None])**2)

# calculate the inverse distances, use a small epsilon to avoid infinities:
idist = 1. / (dist + 1e-12)**pow

# calculate the weighted average for each column with idist as the weight function
unknown_values = np.sum(known_values[None,:] * idist, axis=1) / np.sum(idist, axis=1)
``````

When this example is plotted (blotch area reflects value, red are known points, blue unknown points), we get:

Top