When I'm teaching my courses about scientific Python, students often have a hard time understanding the concept of meshed grids in NumPy and what's the difference between the various options. As this appears to be a quite common problem and as the concept of meshed grids is vital for many scientific applications, I'm writing this post to settle these problems once and for all. ;-)

The Problem

Consider the following code, where we define a one-dimensional equidistant grid x:

import numpy as np

x = np.linspace(-1, 1, 100)

One of the beauties of NumPy is that you can easily apply functions to arrays. For example, we could calculate the sine function on that very grid by

f = np.sin(x)

Both x and f are one-dimensional arrays of length 100, as we would guess:

x.shape

(100,)

f.shape

(100,)

This is all intuitive. But what if we are working in multiple dimensions? Suppose we want to use a two dimensional grid along the \(x\) and \(y\) axes. How do we do that? Specifically, how can we calculate the function

$$ f(x, y) = \sin x \cos y $$

on a two-dimensional grid?

The Solution

First, we define our one-dimensional coordinate values like in the one-dimensional case.

nx = ny = 100
x = np.linspace(-1, 1, nx)
y = np.linspace(-1, 1, ny)

Technically, you could write

python
np.sin(x)*np.cos(y)

but this would only return a one-dimensional array. That's clear because NumPy cannot know that these one-dimensional arrays have the meaning of coordinate values and that the two axes are perpendicular to each other. This is where meshed grids come in.

In a two-dimensional space, a meshgrid consists of two arrays: one representing the x-coordinates and the other representing the y-coordinates. These arrays are created by taking all possible combinations of the x and y values within specified ranges. Each combination of x and y values represents a point on the grid. This creates a rectangular grid-like structure where each point has a unique x and y coordinate.

For example, suppose we want to create a meshgrid for x values ranging from 0 to 2 and y values ranging from 0 to 3. The resulting meshgrid would have two-dimensional X and Y arrays as follows:

X = [[0, 1, 2],
     [0, 1, 2],
     [0, 1, 2],
     [0, 1, 2]]

Y = [[0, 0, 0],
     [1, 1, 1],
     [2, 2, 2],
     [3, 3, 3]]

I usually write the 1-D coordinate arrays in lower case letters, like \(x\), and the meshed grid version in upper case, like \(X\).

Here, each element in the X array corresponds to the x-coordinate of a point on the two-dimensional grid, and each element in the y array corresponds to the y-coordinate of the same point.

The meshgrid function in NumPy simplifies the process of creating meshgrids. It takes as input one or more 1-dimensional arrays representing the coordinates along each dimension and returns a set of arrays representing the meshgrid.

x = [0, 1, 2]
y = [0, 1, 2, 3]

X, Y = np.meshgrid(x, y)

X

array([[0, 1, 2],
       [0, 1, 2],
       [0, 1, 2],
       [0, 1, 2]])

Y

array([[0, 0, 0],
       [1, 1, 1],
       [2, 2, 2],
       [3, 3, 3]])

The output arrays have the same shape, which is the shape of the input arrays broadcasted together:

X.shape

(4, 3)

Y.shape

(4, 3)

With the meshed grid defined, we can now calculate our function on the two-dimensional grid:

f = np.sin(X) * np.cos(Y)

f.shape

(4, 3)

f

array([[ 0.        ,  0.84147098,  0.90929743],
       [ 0.        ,  0.45464871,  0.4912955 ],
       [-0.        , -0.35017549, -0.37840125],
       [-0.        , -0.83304996, -0.90019763]])

Note carefully that the first axis (axis 0) has length 4 and the second axis (axis 1) has length 3. If you compare that with the definition of our one-dimensional coordinate arrays x and y, you notice that obviously meshgrid has put all the \(y\) values into axis 0 of the resulting meshed grid and the \(x\) values into axis 1.

So when you want to the function value for the point \((x, y) = (1, 3)=\) (x[1], y[3]), then we would have to use the syntax

python
f[3, 1]

Depending on your use case, this can be a little unintuitive because we often think of the \(x\)-axis as the first axis and the \(y\)-axis as the second. This strange ordering is the convention in image processing, however, where you think of the first axis as downwards and the second axis as sidewards.

In our case, we would prefer the natural ordering we are used from mathematical functions, instead. We can easily change the behavior of meshgrid by giving it the following indexing='ij' option:

X, Y = np.meshgrid(x, y, indexing='ij')

X

array([[0, 0, 0, 0],
       [1, 1, 1, 1],
       [2, 2, 2, 2]])

Y

array([[0, 1, 2, 3],
       [0, 1, 2, 3],
       [0, 1, 2, 3]])

Now the usage of our grid is much more intuitive for our use case. We again define

f = np.sin(X) * np.cos(Y)

and now we can get the function value for the second point (index 1) on the x axis and the fourth point (index 3) on the y axis intuitively by calling

f[1, 3]

-0.833049961066805

Conclusion

The meshgrid function in NumPy is a practical tool for creating coordinate systems in multiple dimensions. It allows you to generate meshgrid arrays that represent points on a grid, with each point having unique coordinates. By providing one or more 1-dimensional arrays representing the coordinates along each dimension, meshgrid simplifies the process of creating meshgrids.

The default indexing option creates output arrays where the x-coordinates have the same shape as the y-values input array, while the y-coordinates have the same shape as the x-values input array. This convention is suitable applications like image processing.

Alternatively, the 'ij' indexing option can be used to generate output arrays where the x-coordinates have the same shape as the x-values input array, and the y-coordinates have the same shape as the y-values input array. This convention may be preferred in many mathematical contexts.