Plotting contour map of interpolated function: unmatching results for different sections of data

Question

I've determined the value of z = f(x, y) for a grid of points in the plane:

import numpy as np

n_theta, n_phi = 20, 30
thetas = np.linspace(0, 1, n_theta) * np.pi
phis = np.linspace(0, 1, n_phi) * (2 * np.pi)

res = f(thetas, phis)
# 'res' is an np.array with res.shape=(600,3), where each row is [x, y, z]

Given a meshgrid of the x's and y's, I interpolate a function using my previous data for z, and plot a contour map of the result:

from scipy.interpolate import griddata
import matplotlib.pyplot as plt

X, Y = np.meshgrid(res[:,0], res[:,1])
Z = griddata(res[:,0:2], res[:,2], (X, Y), method='cubic')    # Interpolate function.

fig = plt.figure()
contour = plt.contourf(X, Y, Z, levels=100)

This works fine, and the result (example data linked below) is something like this (which isn't actually the shape I was looking for in my particular simulation):

However, when I run the same script with values of phi (horizontal axis) only for half of the original interval (or, equivalently, take only this half of my data), I get a contour like this:

(For the sake of completeness, if I then do it for the second half of the interval and manually put both figures side by side, the result is like the following, which is more like something I would expect to get from the full data):

(If needed, this is the example data to generate the plots, followed by a snippet one could use):

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

# Full axis:

res = np.load("stackoverflow-example-data")
X, Y = np.meshgrid(res[:,0], res[:,1])
Z = griddata(res[:,0:2], res[:,2], (X, Y), method='cubic')    # Interpolate function.
fig = plt.figure()
contour = plt.contourf(X, Y, Z, levels=100)

# Half axis:

half = np.array([res[i,:] for i in range(len(res)) if res[i,0] <= np.pi])
X, Y = np.meshgrid(res[:,0], res[:,1])
Z = griddata(res[:,0:2], res[:,2], (X, Y), method='cubic')    # Interpolate function.
fig = plt.figure()
contour = plt.contourf(X, Y, Z, levels=100)

So what am I missing here?

1. Is this the correct procedure to interpolate and plot the contour as I want?
2. Shouldn't I expect to get similar results for the full axis and only half of it? If not, why?

Answer 1

Your data are already on a grid. Creating a new grid based on that grid generates something very messy.

Now the array res contains 600 x,y,z values. Using reshape(20, 30) tells numpy that these 600 entries are in reality 20 rows of 30 columns. The 'purest' form to show the data, would be a scatterplot, showing just the data.

With imshow the data can be shown as an image. Using interpolation='nearest' would 'blow up' each z-value to fill a rectangular pixel. Using interpolation='bicubic' would smoothly smear out these pixels.

contourf creates contours with equally-valued z. Depending on the data, having many levels (such as 100) would help to get a smooth image, but it wouldn't be valuable compared to just displaying a smoothed image. Having a limited number of levels (e.g. 20) could help to show a general underlying shape.

Here is some code to compare and experiment with the different approaches:

import numpy as np
import matplotlib.pyplot as plt

# Full axis:

res = np.load("stackoverflow-example-data.npy")
X = res[:, 0].reshape(20, 30)
Y = res[:, 1].reshape(20, 30)
Z = res[:, 2].reshape(20, 30)
fig, axes = plt.subplots(ncols=3, figsize=(12, 4))
axes[0].scatter(X, Y, c=Z, s=10)
axes[1].imshow(Z, origin='lower', interpolation='bicubic', aspect='auto',
               extent=[X.min(), X.max(), Y.min(), Y.max()])
contour = axes[2].contourf(X, Y, Z, levels=20)
axes[0].set_title('scatter plot')
axes[1].set_title('imshow, bicubic interpolation')
axes[2].set_title('contour plot, 20 levels')

plt.show()

Experimenting with different colormaps, can accentuate different properties of the function you are showing. You could also transform the Z-value (such as np.log(Z) or np.exp(Z)) so different regions get more or less detail.

for ax, cmap in zip(axes, ('viridis', 'inferno_r', 'Purples')):
    ax.imshow(Z, origin='lower', interpolation='bicubic', aspect='auto', cmap=cmap,
               extent=[X.min(), X.max(), Y.min(), Y.max()])
    ax.set_title(cmap)