Plotting interpolated values using LinearNDInterpolator (Python)

Question

I am using the LinearNDInterpolator on some (x, y, z) data, using the following script. However, I cannot figure out how to go from the interpolated data to plotting/showing the interpolation in heatmap form? Am I missing something like setting up a meshgrid based on the min and max of x and y? Any help or an example would be great!

import numpy as np
import scipy.interpolate

x = np.array([-4386795.73911443, -1239996.25110694, -3974316.43669208,
               1560260.49911342,  4977361.53694849, -1996458.01768192,
               5888021.46423068,  2969439.36068243,   562498.56468588,
               4940040.00457585])

y = np.array([ -572081.11495993, -5663387.07621326,  3841976.34982795,
               3761230.61316845,  -942281.80271223,  5414546.28275767,
               1320445.40098735, -4234503.89305636,  4621185.12249923,
               1172328.8107458 ])

z = np.array([ 4579159.6898615 ,  2649940.2481702 ,  3171358.81564312,
               4892740.54647532,  3862475.79651847,  2707177.605241  ,
               2059175.83411223,  3720138.47529587,  4345385.04025412,
               3847493.83999694])

# Create coordinate pairs
cartcoord = zip(x, y)

# Interpolate
interp = scipy.interpolate.LinearNDInterpolator(cartcoord, z)

Edit: Based on @Spinor's solution, and using Python 2.7, the following code gives me what I'm looking for (approach 1). Is there a way to increase my density of the interpolated points?

The dataset yields the following plot:

Needless to say, I did not expect the results to be circular, since the (lat,lon) coordinates are taken from an equirectrangular projection map. On further investigation, I think this is simply mapped on a different projection.

Answer 1

I will assume that you are trying to interpolate values of z.

Now, what happens when you do call interpolation function? It creates the entire landscape of the inputs (x and y) and the outputs (z). In the code above, you didn't really ask for its value at any point. To use this function, you need to specify the inputs and it will give you the interpolated output.

You used the function scipy.interpolate.LinearNDInterpolator which is constructed by triangulating the input data and on each triangle performing linear barycentric interpolation. Depending on your inputs, there are likely to be regions where this breaks down and you get Nan. For instance, try this in your code

print interp(-4386790, 3720137)

This is within the limits of the min-max of your x and y. We could set the Nan to zero via the fill_value argument if that is acceptable to you.

Read up on the docs. Often people might find the following function acceptable as well, scipy.interpolate.interp2d. It uses spline interpolation instead. In the code below, I've implemented both functions (the former with nan values set to 0) and plotted them on a heatmap.

As for the heatmap, it is as you've suspected. You have to create a grid of values. Below is my the output graphs for LinearNDInterpolator with nan set to zero and interp2d as well as the codes.

Using LinearNDInterpolator(cartcoord, z, fill_value=0)

Using interp2d(x, y, z)

P.S. I am using Python3. If you run into issues in Python2, remove the list from cartcoord = list(zip(x, y)).

import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate

x = np.array([-4386795.73911443, -1239996.25110694, -3974316.43669208,
               1560260.49911342,  4977361.53694849, -1996458.01768192,
               5888021.46423068,  2969439.36068243,   562498.56468588,
               4940040.00457585])

y = np.array([ -572081.11495993, -5663387.07621326,  3841976.34982795,
               3761230.61316845,  -942281.80271223,  5414546.28275767,
               1320445.40098735, -4234503.89305636,  4621185.12249923,
               1172328.8107458 ])

z = np.array([ 4579159.6898615 ,  2649940.2481702 ,  3171358.81564312,
               4892740.54647532,  3862475.79651847,  2707177.605241  ,
               2059175.83411223,  3720138.47529587,  4345385.04025412,
               3847493.83999694])

# Create coordinate pairs
cartcoord = list(zip(x, y))


X = np.linspace(min(x), max(x))
Y = np.linspace(min(y), max(y))
X, Y = np.meshgrid(X, Y)

# Approach 1
interp = scipy.interpolate.LinearNDInterpolator(cartcoord, z, fill_value=0)
Z0 = interp(X, Y)
plt.figure()
plt.pcolormesh(X, Y, Z0)
plt.colorbar() # Color Bar
plt.show()

# Approach 2
func = scipy.interpolate.interp2d(x, y, z)
Z = func(X[0, :], Y[:, 0])
plt.figure()
plt.pcolormesh(X, Y, Z)
plt.colorbar() # Color Bar
plt.show()