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: enter image description here

   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()

What's happening is due to your z-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 do some sanity checks and "whiten" (a.k.a. rescale) the input data, but it doesn't. You should be able to fix your problems by doing a simple linear rescale of `ao` before interpolating and then doing the inverse to the result. — Joe Kington, Jul 10 '13 at 18:38
My guess is that your 'vmin' and 'vmax' values are wrong, I tried it by setting them to 'ao.min()' and 'ao.max()' and it came out like this: http://i.imgur.com/8Om1YdG.png — user1767344, Jul 10 '13 at 18:46
Well, i made a mistake in my earlier comment. It was supposed to be `ai.min()` and `ai.max()` it matters. Or you can simply skip them and matplotlib will fix it automatically. — user1767344, Jul 10 '13 at 18:54

Joe Kington · Accepted Answer · 2013-07-10T19:06:20.607

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()

enter image description here

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.)

score 2 · Answer 2 · answered Aug 19 '13 at 21:13

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()

Matplotlib 2D graph with interpolation

2 Answers2