I need to resample 2D-data to a regular grid.
This is what my code looks like:
import matplotlib.mlab as ml import numpy as np y = np.zeros((512,115)) x = np.zeros((512,115)) # Just random data for this test: data = np.random.randn(512,115) # filling the grid coordinates: for i in range(512): y[i,:]=np.arange(380,380+4*115,4) for i in range(115): x[:,i] = np.linspace(-8,8,512) y[:,i] -= np.linspace(-0.1,0.2,512) # Defining the regular grid y_i = np.arange(380,380+4*115,4) x_i = np.linspace(-8,8,512) resampled_data = ml.griddata(x,y,data,x_i,y_i)
(512,115) is the shape of the 2D data, and I already installed mpl_toolkits.natgrid.
My issue is that I get back a masked array, where most of the entries are nan, instead of an array that is mostly composed of regular entries and just nan at the borders.
Could someone point me to what I am doing wrong?
Comparing your code example to your question's title, I think you're a bit confused...
In your example code, you're creating regularly gridded random data and then resampling it onto another regular grid. You don't have irregular data anywhere in your example...
(Also, the code doesn't run as-is, and you should look into
meshgrid rather than looping through to generate your x & y grids.)
If you're wanting to re-sample an already regularly-sampled grid, as you do in your example, there are more efficient methods than griddata or anything I'm about to describe below. (
scipy.ndimage.map_coordinates would be well suited to your problem, it that case.)
Based on your question, however, it sounds like you have irregularly spaced data that you want to interpolate onto a regular grid.
In that case, you might have some points like this:
import numpy as np import matplotlib.mlab as mlab import matplotlib.pyplot as plt # Bounds and number of the randomly generated data points ndata = 20 xmin, xmax = -8, 8 ymin, ymax = 380, 2428 # Generate random data x = np.random.randint(xmin, xmax, ndata) y = np.random.randint(ymin, ymax, ndata) z = np.random.random(ndata) # Plot the random data points plt.scatter(x,y,c=z) plt.axis([xmin, xmax, ymin, ymax]) plt.colorbar() plt.show()
You can then interpolate the data as you were doing before... (Continued from code snippet above...)
# Size of regular grid ny, nx = 512, 115 # Generate a regular grid to interpolate the data. xi = np.linspace(xmin, xmax, nx) yi = np.linspace(ymin, ymax, ny) xi, yi = np.meshgrid(xi, yi) # Interpolate using delaunay triangularization zi = mlab.griddata(x,y,z,xi,yi) # Plot the results plt.figure() plt.pcolormesh(xi,yi,zi) plt.scatter(x,y,c=z) plt.colorbar() plt.axis([xmin, xmax, ymin, ymax]) plt.show()
However, you'll notice that you're getting lots of artifacts in the grid. This is due to the fact that your x coordinates range from -8 to 8, while y coordinates range from ~300 to ~2500. The interpolation algorithm is trying to make things isotropic, while you may want a highly anisotropic interpolation (so that it appears isotropic when the grid is plotted).
To correct for this, you need to create a new coordinate system to do your interpolation in. There is no one right way to do this. What I'm using below will work, but the "best" way depends heavily on what your data actually represents.
(In other words, use what you know about the system that your data is measuring to decide how to do it. This is always true with interpolation! You should not interpolate unless you know what the result should look like, and are familiar enough with the interpolation algorithm to use that a priori information to your advantage!! There are also much more flexible interpolation algorithms than the Delaunay triangulation that griddata uses by default, as well, but it's fine for a simple example...)
At any rate, one way to do this is to rescale the x and y coordinates so that they range over roughly the same magnitudes. In this case. we'll rescale them from 0 to 1... (forgive the spaghetti string code... I'm just intending this to be an example...)
# (Continued from examples above...) # Normalize coordinate system def normalize_x(data): data = data.astype(np.float) return (data - xmin) / (xmax - xmin) def normalize_y(data): data = data.astype(np.float) return (data - ymin) / (ymax - ymin) x_new, xi_new = normalize_x(x), normalize_x(xi) y_new, yi_new = normalize_y(y), normalize_y(yi) # Interpolate using delaunay triangularization zi = mlab.griddata(x_new, y_new, z, xi_new, yi_new) # Plot the results plt.figure() plt.pcolormesh(xi,yi,zi) plt.scatter(x,y,c=z) plt.colorbar() plt.axis([xmin, xmax, ymin, ymax]) plt.show()
Hope that helps, at any rate... Sorry for the length of the answer!