Resampling irregularly spaced data to a regular grid in Python

Question

Hi!

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?

Thanks!

Answer 1

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.interpolate.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.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!