问题描述:

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

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: