Axes Method - plot_surface¶

The plot_surface method creates a 3D surface plot from grid data. It requires 2D arrays for X, Y coordinates and corresponding Z values.

Understanding meshgrid¶

Before using plot_surface, understand how np.meshgrid creates coordinate grids.

1. Basic meshgrid¶

import numpy as np

x = np.array([0, 1, 2])      # (3,)
y = np.array([-2, -1, 0, 1, 2])  # (5,)
X, Y = np.meshgrid(x, y)
print(X)  # (5, 3)
print(Y)  # (5, 3)

Output:

X:
[[0 1 2]
 [0 1 2]
 [0 1 2]
 [0 1 2]
 [0 1 2]]

Y:
[[-2 -2 -2]
 [-1 -1 -1]
 [ 0  0  0]
 [ 1  1  1]
 [ 2  2  2]]

2. meshgrid for Surface Plots¶

import numpy as np

n = 40
x = np.linspace(-3.0, 3.0, n)  # (40,)
y = np.linspace(-3.0, 3.0, n)  # (40,)
X, Y = np.meshgrid(x, y)       # (40, 40), (40, 40)

Basic Usage¶

Create 3D surface plots with the projection='3d' subplot keyword.

1. Simple Surface¶

import matplotlib.pyplot as plt
import numpy as np

def f(X, Y):
    return X**2 + Y**2

x = np.linspace(-2, 2, 30)
y = np.linspace(-2, 2, 30)
X, Y = np.meshgrid(x, y)
Z = f(X, Y)

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
ax.plot_surface(X, Y, Z)
plt.show()

2. With Colormap¶

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
ax.plot_surface(X, Y, Z, cmap='viridis')
plt.show()

3. With Labels¶

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
ax.plot_surface(X, Y, Z, cmap='viridis')
ax.set_xlabel('X axis')
ax.set_ylabel('Y axis')
ax.set_zlabel('Z axis')
plt.show()

Standard Normal PDF¶

Plot the standard bivariate normal probability density function.

1. Basic Normal PDF¶

import matplotlib.pyplot as plt
import numpy as np

def f(X, Y):
    return np.exp(-X**2 / 2 - Y**2 / 2) / (2 * np.pi)

x = np.linspace(-4, 4, 100)
y = np.linspace(-4, 4, 100)
X, Y = np.meshgrid(x, y)
Z = f(X, Y)

fig, ax = plt.subplots(figsize=(8, 6), subplot_kw={'projection': '3d'})
ax.set_title("Standard Normal PDF", fontsize=15)
ax.plot_surface(X, Y, Z, cmap=plt.cm.coolwarm)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x, y)')
plt.show()

2. With Custom Z-ticks¶

fig, ax = plt.subplots(figsize=(5 * 1.61803398875, 5), subplot_kw={'projection': '3d'})
ax.set_title("Standard Normal PDF", fontsize=15)
ax.plot_surface(
    X, Y, Z,
    rstride=2,
    cstride=2,
    cmap=plt.cm.coolwarm,
    linewidth=0.5,
    antialiased=True
)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x, y)')
ax.set_zticks((0.05, 0.10, 0.15))
plt.show()

Multivariate Normal with Correlation¶

Plot bivariate normal distributions with different correlation coefficients.

1. Setup Function¶

import matplotlib.pyplot as plt
import numpy as np
from scipy import stats

def create_bivariate_normal(mu_1, mu_2, sigma_1, sigma_2, rho):
    return stats.multivariate_normal(
        [mu_1, mu_2],
        [[sigma_1, rho * sigma_1 * sigma_2],
         [rho * sigma_1 * sigma_2, sigma_2]]
    )

2. Single Correlation¶

n = 40
mu_1, mu_2 = 0, 0
sigma_1, sigma_2 = 1, 0.5
rho = 0.0

x = np.linspace(-3.0, 3.0, n)
y = np.linspace(-3.0, 3.0, n)
X, Y = np.meshgrid(x, y)
pos = np.empty(X.shape + (2,))
pos[:, :, 0] = X
pos[:, :, 1] = Y

Z = create_bivariate_normal(mu_1, mu_2, sigma_1, sigma_2, rho)

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
ax.plot_surface(X, Y, Z.pdf(pos), cmap='viridis', linewidth=0)
ax.set_xlabel('X axis')
ax.set_ylabel('Y axis')
plt.show()

3. Comparing Correlations¶

n = 40
mu_1, mu_2 = 0, 0
sigma_1, sigma_2 = 1, 0.5
rho1, rho2, rho3 = 0.0, -0.8, 0.8

x = np.linspace(-3.0, 3.0, n)
y = np.linspace(-3.0, 3.0, n)
X, Y = np.meshgrid(x, y)
pos = np.empty(X.shape + (2,))
pos[:, :, 0] = X
pos[:, :, 1] = Y

Z = lambda rho: stats.multivariate_normal(
    [mu_1, mu_2],
    [[sigma_1, rho * sigma_1 * sigma_2],
     [rho * sigma_1 * sigma_2, sigma_2]]
)

fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(15, 5), subplot_kw={'projection': '3d'})

ax0.plot_surface(X, Y, Z(rho1).pdf(pos), cmap='Oranges', linewidth=0)
ax0.set_xlabel('X axis')
ax0.set_ylabel('Y axis')
ax0.set_title(f'ρ = {rho1}')

ax1.plot_surface(X, Y, Z(rho2).pdf(pos), cmap='viridis', linewidth=0)
ax1.set_xlabel('X axis')
ax1.set_ylabel('Y axis')
ax1.set_title(f'ρ = {rho2}')

ax2.plot_surface(X, Y, Z(rho3).pdf(pos), cmap='PuBuGn', linewidth=0)
ax2.set_xlabel('X axis')
ax2.set_ylabel('Y axis')
ax2.set_title(f'ρ = {rho3}')

plt.tight_layout()
plt.show()

Noisy Surface¶

Add random noise to surface data.

1. Noisy Bowl¶

import matplotlib.pyplot as plt
import numpy as np

def f(X, Y, seed=0):
    np.random.seed(seed)
    return X**2 + Y**2 + 0.15 * np.random.normal(0., 1., X.shape)

x = np.linspace(-1.5, 1.5, 80)
y = np.linspace(-1.5, 1.5, 80)
X, Y = np.meshgrid(x, y)
Z = f(X, Y)

fig, ax = plt.subplots(figsize=(9, 6), subplot_kw={'projection': '3d'})
ax.plot_surface(
    X, Y, Z,
    rstride=2,
    cstride=2,
    cmap=plt.cm.coolwarm,
    linewidth=0.5,
    antialiased=True
)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x, y)')
plt.show()

2. Noise Level Comparison¶

def f_noisy(X, Y, noise_level, seed=0):
    np.random.seed(seed)
    return X**2 + Y**2 + noise_level * np.random.normal(0., 1., X.shape)

x = np.linspace(-1.5, 1.5, 50)
y = np.linspace(-1.5, 1.5, 50)
X, Y = np.meshgrid(x, y)

fig, axes = plt.subplots(1, 3, figsize=(15, 5), subplot_kw={'projection': '3d'})

noise_levels = [0.0, 0.15, 0.5]
titles = ['No Noise', 'Low Noise', 'High Noise']

for ax, noise, title in zip(axes, noise_levels, titles):
    Z = f_noisy(X, Y, noise)
    ax.plot_surface(X, Y, Z, cmap=plt.cm.coolwarm, linewidth=0.5)
    ax.set_title(title)
    ax.set_xlabel('x')
    ax.set_ylabel('y')

plt.tight_layout()
plt.show()

Stride Parameters¶

Control surface mesh density with rstride and cstride.

1. Default Stride¶

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
ax.plot_surface(X, Y, Z, cmap='viridis')
plt.show()

2. Custom Stride¶

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
ax.plot_surface(X, Y, Z, rstride=2, cstride=2, cmap='viridis')
plt.show()

3. Stride Comparison¶

x = np.linspace(-2, 2, 60)
y = np.linspace(-2, 2, 60)
X, Y = np.meshgrid(x, y)
Z = np.sin(np.sqrt(X**2 + Y**2))

fig, axes = plt.subplots(1, 3, figsize=(15, 5), subplot_kw={'projection': '3d'})

strides = [1, 3, 6]

for ax, s in zip(axes, strides):
    ax.plot_surface(X, Y, Z, rstride=s, cstride=s, cmap='viridis')
    ax.set_title(f'rstride={s}, cstride={s}')

plt.tight_layout()
plt.show()

Colormap Options¶

Different colormaps for surface visualization.

1. Sequential Colormaps¶

x = np.linspace(-3, 3, 50)
y = np.linspace(-3, 3, 50)
X, Y = np.meshgrid(x, y)
Z = np.sin(X) * np.cos(Y)

fig, axes = plt.subplots(1, 3, figsize=(15, 5), subplot_kw={'projection': '3d'})

cmaps = ['viridis', 'plasma', 'inferno']

for ax, cmap in zip(axes, cmaps):
    ax.plot_surface(X, Y, Z, cmap=cmap)
    ax.set_title(f"cmap='{cmap}'")

plt.tight_layout()
plt.show()

2. Diverging Colormaps¶

fig, axes = plt.subplots(1, 3, figsize=(15, 5), subplot_kw={'projection': '3d'})

cmaps = ['coolwarm', 'RdBu', 'seismic']

for ax, cmap in zip(axes, cmaps):
    ax.plot_surface(X, Y, Z, cmap=cmap)
    ax.set_title(f"cmap='{cmap}'")

plt.tight_layout()
plt.show()

3. Qualitative Colormaps¶

fig, axes = plt.subplots(1, 3, figsize=(15, 5), subplot_kw={'projection': '3d'})

cmaps = ['Oranges', 'PuBuGn', 'YlOrRd']

for ax, cmap in zip(axes, cmaps):
    ax.plot_surface(X, Y, Z, cmap=cmap)
    ax.set_title(f"cmap='{cmap}'")

plt.tight_layout()
plt.show()

Linewidth and Antialiasing¶

Control edge appearance.

1. No Edge Lines¶

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
ax.plot_surface(X, Y, Z, cmap='viridis', linewidth=0)
plt.show()

2. Visible Edge Lines¶

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
ax.plot_surface(X, Y, Z, cmap='viridis', linewidth=0.5)
plt.show()

3. Antialiasing Comparison¶

fig, axes = plt.subplots(1, 2, figsize=(12, 5), subplot_kw={'projection': '3d'})

axes[0].plot_surface(X, Y, Z, cmap='viridis', linewidth=0.5, antialiased=False)
axes[0].set_title('antialiased=False')

axes[1].plot_surface(X, Y, Z, cmap='viridis', linewidth=0.5, antialiased=True)
axes[1].set_title('antialiased=True')

plt.tight_layout()
plt.show()

Alpha Transparency¶

Control surface transparency.

1. Opaque Surface¶

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
ax.plot_surface(X, Y, Z, cmap='viridis', alpha=1.0)
plt.show()

2. Semi-transparent Surface¶

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
ax.plot_surface(X, Y, Z, cmap='viridis', alpha=0.7)
plt.show()

3. Alpha Comparison¶

fig, axes = plt.subplots(1, 3, figsize=(15, 5), subplot_kw={'projection': '3d'})

alphas = [1.0, 0.7, 0.4]

for ax, a in zip(axes, alphas):
    ax.plot_surface(X, Y, Z, cmap='viridis', alpha=a)
    ax.set_title(f'alpha={a}')

plt.tight_layout()
plt.show()

Common Mathematical Surfaces¶

1. Paraboloid¶

x = np.linspace(-2, 2, 50)
y = np.linspace(-2, 2, 50)
X, Y = np.meshgrid(x, y)
Z = X**2 + Y**2

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
ax.plot_surface(X, Y, Z, cmap='viridis')
ax.set_title('Paraboloid: $z = x^2 + y^2$')
plt.show()

2. Saddle Point¶

Z = X**2 - Y**2

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
ax.plot_surface(X, Y, Z, cmap='coolwarm')
ax.set_title('Saddle: $z = x^2 - y^2$')
plt.show()

3. Sinusoidal Surface¶

Z = np.sin(np.sqrt(X**2 + Y**2))

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
ax.plot_surface(X, Y, Z, cmap='plasma')
ax.set_title('Sinusoidal: $z = sin(\\sqrt{x^2 + y^2})$')
plt.show()

4. Surface Gallery¶

fig, axes = plt.subplots(2, 2, figsize=(12, 10), subplot_kw={'projection': '3d'})

# Paraboloid
axes[0, 0].plot_surface(X, Y, X**2 + Y**2, cmap='viridis')
axes[0, 0].set_title('Paraboloid')

# Saddle
axes[0, 1].plot_surface(X, Y, X**2 - Y**2, cmap='coolwarm')
axes[0, 1].set_title('Saddle')

# Sinusoidal
axes[1, 0].plot_surface(X, Y, np.sin(X) * np.cos(Y), cmap='plasma')
axes[1, 0].set_title('sin(x)cos(y)')

# Gaussian
axes[1, 1].plot_surface(X, Y, np.exp(-(X**2 + Y**2)), cmap='inferno')
axes[1, 1].set_title('Gaussian')

plt.tight_layout()
plt.show()

Full Customization¶

1. Complete Example¶

import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(-3, 3, 80)
y = np.linspace(-3, 3, 80)
X, Y = np.meshgrid(x, y)
Z = np.exp(-(X**2 + Y**2) / 2) / (2 * np.pi)

fig, ax = plt.subplots(figsize=(10, 8), subplot_kw={'projection': '3d'})
surf = ax.plot_surface(
    X, Y, Z,
    rstride=2,
    cstride=2,
    cmap=plt.cm.coolwarm,
    linewidth=0.3,
    antialiased=True,
    alpha=0.9
)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('Bivariate Normal Distribution', fontsize=14)
fig.colorbar(surf, shrink=0.5, aspect=10)
plt.show()

2. Professional Presentation¶

from scipy import stats

n = 50
x = np.linspace(-3, 3, n)
y = np.linspace(-3, 3, n)
X, Y = np.meshgrid(x, y)
pos = np.dstack((X, Y))

rv = stats.multivariate_normal([0, 0], [[1, 0.5], [0.5, 1]])
Z = rv.pdf(pos)

fig, ax = plt.subplots(figsize=(10, 8), subplot_kw={'projection': '3d'})
surf = ax.plot_surface(
    X, Y, Z,
    cmap='viridis',
    linewidth=0,
    antialiased=True
)
ax.set_xlabel('X', fontsize=12)
ax.set_ylabel('Y', fontsize=12)
ax.set_zlabel('Density', fontsize=12)
ax.set_title('Bivariate Normal (ρ = 0.5)', fontsize=14, fontweight='bold')
fig.colorbar(surf, shrink=0.6, aspect=15, label='Probability Density')
plt.tight_layout()
plt.show()