use of matlab meshgrid

Question

Sorry i have to keep coming back to the problem gridding irregularly spaced data. I do not seem to see any clear responses to questions of how to grid data to a regular grid and the software documentation to me is good for those who already know. I have x, y, z data on 29 points, with a header "Lon Lat Z". to plot contours with this data here is what I do:

1. After reading in the data, make a 300 by 300 point regular grid onto which to interpolate

   numcols, numrows = 300, 300
   xi = np.linspace(data.Lon.min(), data.Lon.max(), numcols)
   yi = np.linspace(data.Lat.min(), data.Lat.max(), numcols)
   xi, yi = np.meshgrid(xi, yi)
   

   Print xi and print yi at this point gives me x and y according to my data, interpolated over 300x300 points.

2. Interpolate the data over the grid created above

   x, y, z = data.Lon.values, data.Lat.values, data.Z.values
   zi = griddata(x, y, z, xi, yi)
   

   At this point if I do print zi I get

   [[-- -- -- ..., -- -- --]
   [-- -- -- ..., -- -- --]
   [-- -- -- ..., -- -- --]
   ...,
   [-- -- -- ..., -- -- --]
   [-- -- -- ..., -- -- --]
   [-- -- -- ..., -- -- --]]

I was expecting to see values of interpolated z. I also have a map object defined to be overlaid by contours. The plotting function gives me separate figures, for the basemap and for the contours, with the correct contour values. My question is why am I getting blank values for contours and how come they are plotted correctly? For completeness here is my plotting function

fig=plt.figure(figsize=(8,4.5))
im = plt.contourf(xi, yi, zi)
plt.show()

Two plots come up (a base map and contours side by side)

Any help please.

Answer 1

The following should work:

numcols, numrows = 300, 300
xi = np.linspace(data.Lon.min(), data.Lon.max(), numrows)
yi = np.linspace(data.Lat.min(), data.Lat.max(), numcols)
xi, yi = np.meshgrid(xi, yi)

x, y, z = data.Lon.values, data.Lat.values, data.Z.values
points = np.vstack((x,y)).T
values = z
wanted = (xi, yi)
zi = griddata(points, values, wanted)

So that last line is how griddata works (assuming you use scipy.interpolate.griddata?) The problem you have is that you seem to give griddata five arguments, while if i look at http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html#scipy.interpolate.griddata It says the following:

scipy.interpolate.griddata(points, values, xi, method='linear', fill_value=nan)

So in your case, giving five arguments is where it goes wrong I guess (confirm if it works if you do it like this, since I don't have your data so I can't test if it gives the correct result).

So in your case, your points where the values are known are the x, while the values which are known at that point are the y-values, and the points where you want to know them are at the z-values. No idea how the method='linear' copes with your argument, and the fill_value you give in is bad too, so you should just give the right inputs (which I think are correct the way I formulated them), and then it should work right.

edit: read in your data as a txt, and wrote the following code. Can you run it to see if that is the result you wanted?

import numpy as np
from scipy.interpolate import griddata
class d():
    def __init__(self):
        A0 = open("test.txt","rb") # i just copypasted your data into a txt (without first row), and reading it in in this class, so that the names are the same as yours
        A1 = A0.readlines()
        A = np.zeros((len(A1),3))
        for i, l in enumerate(A1):
            li = l.split()
            A[i,0] = float(li[0])
            A[i,1] = float(li[1])
            A[i,2] = float(li[2])
        self.Lon = A[:,0]
        self.Lat = A[:,1]
        self.Z = A[:,2]

data = d()
numcols, numrows = 30, 30
xi = np.linspace(data.Lon.min(), data.Lon.max(), numrows)
yi = np.linspace(data.Lat.min(), data.Lat.max(), numcols)
xi, yi = np.meshgrid(xi, yi)

x, y, z = data.Lon, data.Lat, data.Z
points = np.vstack((x,y)).T
values = z
wanted = (xi, yi)
zi = griddata(points, values, wanted)
import pylab as plt
fig = plt.figure(0, figsize=(8,4.5))
im = plt.contourf(xi, yi, zi)
plt.colorbar()
fig2 = plt.figure(1, figsize=(8,4.5))
im = plt.scatter(xi, yi, c= zi)
plt.colorbar()
plt.show()