NumPy - np.meshgrid() with Examples
import numpy as np
Key Insights
np.meshgrid()creates coordinate matrices from coordinate vectors, essential for vectorized function evaluation over multi-dimensional grids without explicit loops- Understanding the
indexingparameter (‘xy’ vs ‘ij’) prevents coordinate system confusion when working with mathematical functions versus matrix operations - Meshgrid enables efficient computation of 3D surfaces, contour plots, and numerical methods like gradient descent by broadcasting operations across entire grids
What np.meshgrid() Does
np.meshgrid() transforms one-dimensional coordinate arrays into multi-dimensional coordinate grids. Instead of writing nested loops to evaluate functions at every point combination, meshgrid creates matrices where each element represents coordinates in N-dimensional space.
import numpy as np
x = np.array([1, 2, 3])
y = np.array([4, 5])
X, Y = np.meshgrid(x, y)
print("X:\n", X)
print("\nY:\n", Y)
Output:
X:
[[1 2 3]
[1 2 3]]
Y:
[[4 4 4]
[5 5 5]]
Each point (X[i,j], Y[i,j]) represents a coordinate pair in the grid. The first array repeats x-values across rows, while the second repeats y-values down columns.
Basic 2D Grid Construction
Creating evaluation grids for mathematical functions demonstrates the core use case:
import numpy as np
import matplotlib.pyplot as plt
# Define coordinate ranges
x = np.linspace(-5, 5, 50)
y = np.linspace(-5, 5, 50)
# Create meshgrid
X, Y = np.meshgrid(x, y)
# Evaluate function: z = x^2 + y^2
Z = X**2 + Y**2
print(f"Grid shape: {X.shape}")
print(f"Total points: {X.size}")
print(f"First few X coordinates:\n{X[:3, :3]}")
print(f"First few Y coordinates:\n{Y[:3, :3]}")
This creates a 50×50 grid (2,500 points) without loops. Each operation on X and Y applies element-wise across the entire grid through NumPy broadcasting.
Indexing Parameter: ‘xy’ vs ‘ij’
The indexing parameter controls coordinate system orientation. Default 'xy' uses Cartesian indexing (standard for plotting), while 'ij' uses matrix indexing.
import numpy as np
x = np.array([1, 2, 3])
y = np.array([4, 5, 6])
# Cartesian indexing (default)
X_cart, Y_cart = np.meshgrid(x, y, indexing='xy')
# Matrix indexing
X_mat, Y_mat = np.meshgrid(x, y, indexing='ij')
print("Cartesian indexing (xy):")
print("X shape:", X_cart.shape) # (3, 3) - rows vary with y
print("X:\n", X_cart)
print("Y:\n", Y_cart)
print("\nMatrix indexing (ij):")
print("X shape:", X_mat.shape) # (3, 3) - rows vary with x
print("X:\n", X_mat)
print("Y:\n", Y_mat)
Output shows the transpose relationship:
Cartesian indexing (xy):
X shape: (3, 3)
X:
[[1 2 3]
[1 2 3]
[1 2 3]]
Y:
[[4 4 4]
[5 5 5]
[6 6 6]]
Matrix indexing (ij):
X shape: (3, 3)
X:
[[1 1 1]
[2 2 2]
[3 3 3]]
Y:
[[4 5 6]
[4 5 6]
[4 5 6]]
Use 'xy' for plotting and mathematical functions. Use 'ij' when working with image processing or matrix operations where row/column semantics matter.
3D Surface Plotting
Meshgrid excels at generating data for 3D visualizations:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Create grid
x = np.linspace(-3, 3, 100)
y = np.linspace(-3, 3, 100)
X, Y = np.meshgrid(x, y)
# Gaussian surface
Z = np.exp(-(X**2 + Y**2) / 2)
# Plot
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(X, Y, Z, cmap='viridis')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.colorbar(surf)
plt.show()
The same pattern works for complex functions:
# Saddle surface
Z_saddle = X**2 - Y**2
# Ripple effect
Z_ripple = np.sin(np.sqrt(X**2 + Y**2))
# Combined function
Z_complex = np.sin(X) * np.cos(Y) * np.exp(-(X**2 + Y**2) / 10)
Contour Plots and Level Sets
Meshgrid enables contour plotting for visualizing function behavior:
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-2, 2, 200)
y = np.linspace(-2, 2, 200)
X, Y = np.meshgrid(x, y)
# Rosenbrock function (optimization test function)
a, b = 1, 100
Z = (a - X)**2 + b * (Y - X**2)**2
# Create contour plot
plt.figure(figsize=(10, 8))
contours = plt.contour(X, Y, Z, levels=20, cmap='RdYlBu_r')
plt.clabel(contours, inline=True, fontsize=8)
plt.colorbar(contours)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Rosenbrock Function Contours')
plt.grid(True, alpha=0.3)
plt.show()
This visualization helps understand optimization landscapes, gradient directions, and function minima.
Multi-Dimensional Grids
Meshgrid handles arbitrary dimensions through the sparse parameter:
import numpy as np
# 3D grid
x = np.linspace(0, 1, 3)
y = np.linspace(0, 1, 4)
z = np.linspace(0, 1, 5)
X, Y, Z = np.meshgrid(x, y, z, indexing='ij')
print(f"3D grid shapes: X={X.shape}, Y={Y.shape}, Z={Z.shape}")
# Evaluate 3D function
result = X**2 + Y**2 + Z**2
print(f"Result shape: {result.shape}")
# Sparse grid (memory efficient)
X_sparse, Y_sparse, Z_sparse = np.meshgrid(x, y, z, indexing='ij', sparse=True)
print(f"\nSparse shapes: X={X_sparse.shape}, Y={Y_sparse.shape}, Z={Z_sparse.shape}")
# Same computation works due to broadcasting
result_sparse = X_sparse**2 + Y_sparse**2 + Z_sparse**2
print(f"Sparse result shape: {result_sparse.shape}")
print(f"Results equal: {np.allclose(result, result_sparse)}")
Sparse grids store only one-dimensional arrays, reducing memory from O(n^d) to O(n*d) for d dimensions with n points per axis.
Vectorized Distance Calculations
Computing distances from every grid point to specific locations:
import numpy as np
# Create grid
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x, y)
# Multiple point sources
sources = np.array([[1, 1], [-2, 2], [0, -3]])
# Calculate distance fields
distance_fields = []
for source in sources:
dist = np.sqrt((X - source[0])**2 + (Y - source[1])**2)
distance_fields.append(dist)
# Minimum distance to any source
min_distance = np.minimum.reduce(distance_fields)
# Find nearest source for each grid point
nearest_source = np.argmin([df for df in distance_fields], axis=0)
print(f"Distance field shape: {min_distance.shape}")
print(f"Min distance: {min_distance.min():.4f}")
print(f"Max distance: {min_distance.max():.4f}")
This pattern applies to Voronoi diagrams, nearest neighbor searches, and field simulations.
Numerical Gradient Computation
Meshgrid facilitates gradient calculations for optimization:
import numpy as np
def function(x, y):
return x**2 + 2*y**2 - 2*x*y + 2*x - 4*y
# Create fine grid
x = np.linspace(-5, 5, 50)
y = np.linspace(-5, 5, 50)
X, Y = np.meshgrid(x, y)
# Evaluate function
Z = function(X, Y)
# Compute gradients using numpy
grad_y, grad_x = np.gradient(Z, y, x)
# Find minimum (where gradient is near zero)
grad_magnitude = np.sqrt(grad_x**2 + grad_y**2)
min_idx = np.unravel_index(np.argmin(grad_magnitude), grad_magnitude.shape)
print(f"Approximate minimum at: x={X[min_idx]:.3f}, y={Y[min_idx]:.3f}")
print(f"Function value: {Z[min_idx]:.3f}")
print(f"Gradient magnitude: {grad_magnitude[min_idx]:.6f}")
This approach visualizes gradient descent paths and identifies critical points in optimization problems.
Performance Considerations
Meshgrid trades memory for speed. For large grids, consider alternatives:
import numpy as np
import time
# Large grid comparison
n = 5000
# Method 1: meshgrid (memory intensive)
start = time.time()
x = np.linspace(0, 1, n)
y = np.linspace(0, 1, n)
X, Y = np.meshgrid(x, y)
Z1 = X**2 + Y**2
time_meshgrid = time.time() - start
# Method 2: sparse meshgrid
start = time.time()
X_sparse, Y_sparse = np.meshgrid(x, y, sparse=True)
Z2 = X_sparse**2 + Y_sparse**2
time_sparse = time.time() - start
# Method 3: broadcasting without meshgrid
start = time.time()
Z3 = x[:, np.newaxis]**2 + y[np.newaxis, :]**2
time_broadcast = time.time() - start
print(f"Meshgrid time: {time_meshgrid:.4f}s")
print(f"Sparse time: {time_sparse:.4f}s")
print(f"Broadcasting time: {time_broadcast:.4f}s")
print(f"All methods equal: {np.allclose(Z1, Z2) and np.allclose(Z1, Z3)}")
For simple operations, direct broadcasting often matches meshgrid performance while using less memory. Use meshgrid when code clarity matters or when passing grids to functions expecting explicit coordinate arrays.