[FIXED] Radially 'sweep out' a 2D contour plot to create a 3D plot (Python, Matplotlib)

Issue

I have three 1D arrays, which represent radius, height, and an intensity measured at that point. I have plotted these to create a 2D contour map. A simple example of the way in which the data is stored is below:

import numpy as np
import matplotlib.pyplot as plt

x = [1,1,1,2,2,2,3,3,3]
y = [1,2,3,1,2,3,1,2,3]
intensity = [5,6,8,9,9,11,15,5,2]

plt.xlabel('Radius')
plt.ylabel('Height')

plt.tricontourf(x,y,intensity)
plt.colorbar(label='Intensity')
plt.show()

(I have had to use plt.tricontourf rather than plt.contour, since the z data is not 2D)

I am looking to create a 3D plot by 'sweeping' the 2D plot through 360 degrees, creating a disk which is azimuthally symmetric. The image below illustrates what I am trying to do...

...with the data interpolated smoothly through the 360 degrees.

There are a couple of similar questions, notably this one, but this does not use three sets of data to create the contours.

Solution

Technically you cannot rotate a 2D plot and get a 3D surface. You can only rotate a 2D curve and get a 3D surface. If this is the case, you could do it as:

import numpy as np
from matplotlib import cm
import matplotlib.pyplot as plt


fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(projection='3d')
N = 100
r = np.linspace(0, 1, N)
z = np.sqrt(1 - r**2)
intensity = np.linspace(0, 1, N).reshape(1, -1)

theta = np.linspace(0, 2*np.pi-1e-3, N)
X = np.outer(np.cos(theta), r)
Y = np.outer(np.sin(theta), r)
Z = np.repeat(z.reshape(1, -1), N, axis = 0)

surf = ax.plot_surface(X, Y, Z, facecolors=cm.jet(np.repeat(intensity, N, axis = 0)))
ax.axes.set_zlim3d(-1, 1)
plt.show()

In the code I rotated a curve to create half a unit sphere and color it according to intensity:

to

If you insist on plotting all the points, I would suggest a 3d scatter plot, I did some linear interpolation to show more points than the original 9:

from scipy.interpolate import interp2d
x = [1,1,1,2,2,2,3,3,3]
y = [1,2,3,1,2,3,1,2,3]
intensity = [5,6,8,9,9,11,15,5,2]
# number of points to interpolate in 3d space
N = 36
# number of points to interpolate in 2d space
N_2d = 10

f = interp2d(x, y, intensity)
# sample along the radius
r = np.linspace(1,3,N_2d)
# sample along z
z = np.linspace(1,3,N_2d)
intensity = f(r, z)
r,z = np.meshgrid(r, z)
theta = np.linspace(0, 2*np.pi, N)
X = np.outer(np.cos(theta), r)
Y = np.outer(np.sin(theta), r)
Z = np.repeat(z.reshape(1, -1), N, axis = 0)

fig = plt.figure(figsize = (10, 6))
ax = fig.add_subplot(projection='3d')
ax.scatter3D(X, Y, Z, c=np.tile(intensity.T, N).T, alpha = 0.5)
plt.show()

Answered By - Z Li