0

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  v7
10  10  1   2   4   3   NA      4
10  15  4   4   3   NA  NA      NA
30  35  NA  NA  1   NA  5       18
5   20  4   NA  4   3   10      NA
15  15  NA  5   4   NA  NA      5
25  10  7   8   7   5   10      NA

My grid point data looks

grid.txt

 X  Y
10  10
10  15
10  20
10  25
10  30
10  35
10  40
15  10
15  15
15  20
15  25
15  30
15  35
15  40
20  10
20  15
20  20
20  25
20  30
20  35
20  40
25  10
25  15
25  20
25  25
25  30
25  35
25  40
30  10
30  15
30  20
30  25
30  30
30  35
30  40
35  10
35  15
35  20
35  25
35  30
35  35
35  40
40  10
40  15
40  20
40  25
40  30
40  35
40  40

The code I am trying is

#! /usr/bin/python

 from math import pow  
 from math import sqrt  
 import numpy as np  

 def pointValue(x,y,power,smoothing,xv,yv,values):  
     nominator=0  
     denominator=0  
     for 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 singularities  
         if(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 zero  
     if denominator > 0:  
         value = nominator/denominator  
     else:  
         value = -9999  
     return value  

 def 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 valuesGrid  


 if __name__ == "__main__":  
     power=1  
     smoothing=20  

     #Creating some data, with each coodinate and the values stored in separated lists  
     xv = [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 value  
     ZI = invDist(xv,yv,values,100,100,power,smoothing)  
     print ZI

1 Answer 1

0

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:

enter image description here

1
  • Thank you for you solution DrV. I will have a look at the other possibilities as well. I am very new to python and trying to learn it hard. My main problems are how to do the loop to iterate though my data (in my example data file, I have more than one column of data to be interpolated), how to do the input output ( because my real data is very big so that I want the code to read both the station data and grid point data as i tried to sohw in the example), and how to ignore the data with null values for the calculation.
    – AiE
    Aug 8, 2014 at 21:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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