Saturday, March 3, 2012

Creating a grid from scattered data using inverse of the distance with python


Attention: I've found some problems in this post. Please read the new one about the same topic.


A usual problem is drawing a map from some scattered data. So it's necessary to set a value for each point in the map from the data in the points already known.
There are different methods used for this purpose, such as kriging or inverse of the distance. The second one is the one used in the example.
The python library scipy has a function called RBF that does that. The code:


import numpy as np
from scipy.interpolate import Rbf

import matplotlib.pyplot as plt

#Creating some data, with each coordinate and the values stored in separated lists
x = [10,60,40,70,10,50,20,70,30,60]
y = [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.0, 100)
XI, YI = np.meshgrid(ti, ti)

#Creating the interpolation function and populating the output matrix value
rbf = Rbf(x, y, values, function='inverse')
ZI = rbf(XI, YI)

# Plotting the result
n = plt.normalize(0.0, 100.0)
plt.subplot(1, 1, 1)
plt.pcolor(XI, YI, ZI)
plt.scatter(x, y, 100, values)
plt.title('RBF interpolation')
plt.xlim(0, 100)
plt.ylim(0, 100)
plt.colorbar()

plt.show() 

When executing the script, the output is an image like this:
The Rbf function  takes as arguments the x and y axes, and a list of the values in the points.  A function to use has to be set as well. In our example, inverse of the distance is used, but there are other options, and even a custom function can be used.
Then an output matrix has to be defined. To do it, the new x and y axes must be defined. In our example, a 100x100 pixel matrix has been defined with the meshgrid function.
Finally, the output matrix is created with the function that was th output of the Rbf function. This ZI matrix is what we can use to draw a map.
In this example, there is only a simple plot, but in the next post a real case will be shown.

1 comment:

  1. How do you do a kriging interpolation? Please help me

    ReplyDelete