Matplotlib 2D graph with interpolation

Question

I'm trying to create a plot to show a temperature gradient over a level surface. I am able to get the colors down but it is not interpolating the way I need it to. This question was helpful in creating the plot. It seems as if the points aren't weighted equally when the plot is produced. This is how the graph currently looks. It should look like the color bar, since to top two points are 23.7 and the bottom two points are 23.4. The code I used is:

   xi, yi = np.linspace(xm, xM, 500), np.linspace(ym, yM, 500)
   xi, yi = np.meshgrid(xi, yi)
   xo=[0, 0, 2, 2]
   yo=[0, 2, 0, 2]
   ao=[23.4, 23.7, 23.4, 23.7]
   rbf = scipy.interpolate.Rbf(xo, yo, ao, function='linear')
   ai = rbf(xi, yi)

   plt.xlabel('X')
   plt.ylabel('Y')

   plt.imshow(ai, vmin=am, vmax=aM, origin='lower', extent=[xm, xM, ym, yM])

   plt.scatter(xo, yo, c=ao)
   plt.colorbar()
   plt.show()

Solution 1

To expand on my earlier comment:

What's happening is due to your a-values having absolute values ~10x greater than the x and y values.

This leads to numerical stability issues when solving for the appropriate weights. Scipy should probably do some sanity checks and "whiten" (a.k.a. rescale) the input data, but it doesn't.

Due to the incorrect weights, the interpolated values are below the mininimum input a values (which shouldn't happen with a linear RBF). This leads to your specified vmin and vmax values "clipping" the grid, as @user1767344 noticed. I'm displaying the "unclipped" version below, but if you specify a similar vmin and vmax, you'll see the same results as your original example.

As an example:

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

def main():
    x = np.array([0, 0, 2, 2])
    y = np.array([0, 2, 0, 2])
    a = np.array([23.4, 23.7, 23.4, 23.7])
    xi, yi = np.mgrid[x.min():x.max():500j, y.min():y.max():500j]

    a_orig = normal_interp(x, y, a, xi, yi)
    a_rescale = rescaled_interp(x, y, a, xi, yi)

    plot(x, y, a, a_orig, 'Not Rescaled')
    plot(x, y, a, a_rescale, 'Rescaled')
    plt.show()

def normal_interp(x, y, a, xi, yi):
    rbf = scipy.interpolate.Rbf(x, y, a)
    ai = rbf(xi, yi)
    return ai

def rescaled_interp(x, y, a, xi, yi):
    a_rescaled = (a - a.min()) / a.ptp()
    ai = normal_interp(x, y, a_rescaled, xi, yi)
    ai = a.ptp() * ai + a.min()
    return ai

def plot(x, y, a, ai, title):
    fig, ax = plt.subplots()

    im = ax.imshow(ai.T, origin='lower',
                   extent=[x.min(), x.max(), y.min(), y.max()])
    ax.scatter(x, y, c=a)

    ax.set(xlabel='X', ylabel='Y', title=title)
    fig.colorbar(im)

main()      

The only difference between the two is just to linearly rescale the input and output a values between 0 and 1:

def normal_interp(x, y, a, xi, yi):
    rbf = scipy.interpolate.Rbf(x, y, a)
    ai = rbf(xi, yi)
    return ai

def rescaled_interp(x, y, a, xi, yi):
    a_rescaled = (a - a.min()) / a.ptp()
    ai = normal_interp(x, y, a_rescaled, xi, yi)
    ai = a.ptp() * ai + a.min()
    return ai

Because it's a 2D interpolation method, the result isn't exactly "like the colorbar". If you want it be, you can just use a 1D interpolation and tile the results. Alternately, this is a case where simple triangulation interpolation methods (e.g. griddata) are a good choice, and should give a result identical to the 1D result. (The disadvantage is that it won't be "smooth" in other cases, but that's always a tradeoff.)

Solution 2

Alternatively, if you have installed matplolib 1.3.0, you can triangulate your data and use a refining / smoothing function : matplotlib.tri.UniformTriRefiner (it should provide smooth contours for any a data and is exact for linear functions of x, y)

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
import matplotlib.cm as cm

def main():
    x = np.array([0, 0, 2, 2])
    y = np.array([0, 2, 0, 2])
    a = np.array([23.4, 23.7, 23.4, 23.7])

    triang = mtri.Triangulation(x, y)
    refiner = mtri.UniformTriRefiner(triang)
    tri_refi, z_test_refi = refiner.refine_field(a, subdiv=4)

    plt.figure()
    plt.gca().set_aspect('equal')
    levels = np.arange(23.4, 23.7, 0.025)
    cmap = cm.get_cmap(name='terrain')
    plt.tricontourf(tri_refi, z_test_refi, levels=levels, cmap=cmap)
    plt.colorbar()

    plt.title("Triangulated")
    plt.show()

main()