You have come across** numpy.meshgrid() **already. You may even have used it. But do you know what it does and how it does it?

If your answer is “I’m not sure”, then you’re not alone. Some people find this function hard to understand. Others understand what it does but not why it’s needed. And some others understand why it’s needed but not how it works.

The first time I came across meshgrid wasn’t in Python or NumPy but in another programming language called MATLAB which I used when I worked as a research scientist.

I was learning programming ‘on the job’ during my PhD studies. So MATLAB’s meshgrid was one of the very first functions I came across in my programming journey as it was essential for my research work at the time. I remember clearly how, for a long time, I used it in my code but didn’t really get it.

When, years later, I came across Python’s `numpy.meshgrid()`

, I was lucky enough to have used its MATLAB counterpart for many years.

In this article, you’ll explore both how `numpy.meshgrid()`

works and when you’ll need it. You’ll see alternatives to using `meshgrid()`

that may be better in some instances.

*If you’re already familiar with the basics of NumPy, you can skip directly to the section ‘Converting The Sine Into A 2D Mathematical Function’.*

## Setting The Scene For `numpy.meshgrid()`

You’ll use NumPy and Matplotlib in this article. You’re likely to have used these packages already if you’re reading this post. If you haven’t, then you’ll need in install them by using `pip`

in the terminal:

$ python -m pip install numpy $ python -m pip install matplotlib

If you’re not familiar with these packages, you can read *Numerical Python for Quantitative Applications Using NumPy* and *Basics of Data Visualisation in Python Using Matplotlib*, two of the chapters in The Python Coding Book.

You can now open a new Python script called `exploring_meshgrid.py`

in your editor and import NumPy and Matplotlib:

# exploring_meshgrid.py import matplotlib.pyplot as plt import numpy as np x = np.linspace(-1, 1, 1000) print(x)

You create an array using `np.linspace()`

. When you run this script, you’ll see this one-dimensional array printed out with values ranging from `-1`

to `1`

. The array has `1,000`

points:

[-1. -0.997998 -0.995996 -0.99399399 -0.99199199 -0.98998999 -0.98798799 -0.98598599 -0.98398398 -0.98198198 -0.97997998 -0.97797798 -0.97597598 -0.97397397 -0.97197197 -0.96996997 -0.96796797 -0.96596597 ... 0.96996997 0.97197197 0.97397397 0.97597598 0.97797798 0.97997998 0.98198198 0.98398398 0.98598599 0.98798799 0.98998999 0.99199199 0.99399399 0.995996 0.997998 1. ]

I truncated the output above. You’ll see the full array when you run the script. NumPy’s `linspace()`

creates a **linear space**. This is a 1D array with equally spaced points.

You can read more about `np.linspace()`

in the very first article I ever wrote for Real Python: `np.linspace()`

: Create Evenly or Non-Evenly Spaced Arrays.

Although we usually prefer to avoid single-letter variable names and use more descriptive names, in this instance, I’m using the variable name `x`

as this mirrors the mathematical convention of using *x* in the Cartesian coordinate system.

## Creating A 1D Mathematical Function

You can use the array `x`

to represent 1D mathematical functions. A quick reminder that a *function* in maths is not the same entity as a *function* in programming. In this article, I use the term *mathematical function* at times but I also use the term *function* by itself at other times when the context will make it clear whether I’m referring to a Python function or a mathematical function.

You’ll plot the following function:

y = \sin(2\pi x/\lambda)

This is a sine wave with wavelength \lambda:

# exploring_meshgrid.py import matplotlib.pyplot as plt import numpy as np wavelength = 0.5 x = np.linspace(-1, 1, 1000) y = np.sin(2 * np.pi * x / wavelength) plt.plot(x, y) plt.show()

When you run this code, you’ll get a sine wave:

## Converting The Sine Into A 2D Mathematical Function

You can write the 1D equation you just plotted using the mathematical functional notation:

f(x) = \sin(2\pi x/\lambda)

This shows more explicitly that this is a function of *x*. The function’s value depends on the position of *x*. You can convert this into a 2D equation by using a function of *x* and *y*:

f(x, y) = \sin(2\pi x/\lambda)

In this example, the right-hand side hasn’t changed. However, the left-hand side shows that this is a 2D equation. In this case, the function only depends on the value of *x* and it is constant in *y*.

You can see what this looks like by plotting using Python. But first, you need to find a way of representing a 2D equation using NumPy.

## Going From 1D To 2D Using `numpy.meshgrid()`

The function *f(x, y)* is a 2D function. Therefore, you’ll need a 2D array to represent its values in Python.

However, the variable `x`

is a 1D array. When you use `x`

as an argument for `np.sin()`

, the result is another 1D array. This is the 1D example you saw earlier.

How can you convince Python that, this time, you’d like to have a 2D function as a result?

You’ll need a version of `x`

that is 2D instead of 1D. What does such a 2D array look like?

Here’s a representation of the 1D variable `x`

showing just `5`

values:

Its 2D counterpart is uppercase `X`

and looks like this:

The values from `0`

to `4`

are repeated for each value along the *y*-axis. I’m showing a square array in this case. However, the *x-* and *y-*dimensions don’t need to be identical.

You can also create an array to represent values along the *y-*axis before creating the 2D versions. I’m choosing different values for this `y`

array for clarity:

This 1D array’s 2D counterpart is `Y`

and looks like this:

The values `12`

to `16`

are now repeated for each value along the *x-*axis.

I’m using NumPy’s `"xy"`

indexing system which is the Cartesian coordinate system. The first index represents columns and changes along the horizontal axis. The second index represents rows and changes along the vertical axis. Note that this Cartesian system is flipped along the horizontal and values increase downwards for *y*.

The alternative is to use the `"ij"`

indexing system which represents standard matrix indexing. The first index represents rows and the second represents columns in this system.

We’ll talk a bit more about the two indexing options later on when moving to higher dimensions.

### Using `numpy.meshgrid()`

How can you create the 2D versions from the 1D vectors? That’s where `numpy.meshgrid()`

enters the stage. You can explore this in the console/REPL first:

>>> import numpy as np >>> x = np.array(range(5)) >>> x array([0, 1, 2, 3, 4]) >>> y = np.array(range(12, 17)) >>> y array([12, 13, 14, 15, 16]) >>> X, Y = np.meshgrid(x, y) >>> X array([[0, 1, 2, 3, 4], [0, 1, 2, 3, 4], [0, 1, 2, 3, 4], [0, 1, 2, 3, 4], [0, 1, 2, 3, 4]]) >>> Y array([[12, 12, 12, 12, 12], [13, 13, 13, 13, 13], [14, 14, 14, 14, 14], [15, 15, 15, 15, 15], [16, 16, 16, 16, 16]])

In this example, you first create the `x`

and `y`

vectors you saw in the diagrams in the previous section. You use `numpy.meshgrid()`

to create the 2D arrays `X`

and `Y`

. These arrays correspond to the previous section’s diagrams for `X`

and `Y`

.

By default, `numpy.meshgrid()`

uses the `"xy"`

indexing system. This will do as I’ll be using `meshgrid()`

to create equations in the Cartesian coordinate system. Choose the indexing system that makes more sense for your application.

There’s one more thing before returning to the 2D version of the sine wave. If you’re a disciplined Python programmer, you may be wondering why I broke with Python convention and used uppercase for the Python variable names `X`

and `Y`

. This is a convention used in MATLAB’s meshgrid and elsewhere.

After so many years of using MATLAB’s meshgrid, I can’t get myself to change this habit. It’s one of the very rare cases where I choose to go against Python’s naming conventions. Many other programmers also stick with uppercase letters to represent the 2D versions of vectors returned by `meshgrid()`

in Python.

### Creating The 2D Sine Wave

You’re now ready to plot this function:

f(x, y) = \sin(2\pi x/\lambda)

All you need to do in the Python code is replace the 1D vector `x`

with its 2D counterpart `X`

. First, you’ll need to call `numpy.meshgrid()`

:

# exploring_meshgrid.py import matplotlib.pyplot as plt import numpy as np wavelength = 0.5 x = y = np.linspace(-1, 1, 1000) X, Y = np.meshgrid(x, y) output = np.sin(2 * np.pi * X / wavelength) plt.imshow(output) plt.show()

This gives the following output, which shows a 2D version of the sine wave:

You can change the colour map to grey and correct the labels on the *x-* and *y-*axes:

# exploring_meshgrid.py import matplotlib.pyplot as plt import numpy as np wavelength = 0.5 x = y = np.linspace(-1, 1, 1000) X, Y = np.meshgrid(x, y) output = np.sin(2 * np.pi * X / wavelength) plt.imshow( output, cmap="gray", extent=[np.min(x), np.max(x), np.min(y), np.max(y)] ) plt.show()

The 2D sine wave now looks like this:

You use `meshgrid()`

to convert the 1D vectors representing the axes into 2D arrays. You can then use those arrays in place of the *x* and *y* variables in the mathematical equation. Since `X`

is a 2D NumPy array, you’ll get a 2D array back when you use `X`

in `np.sin()`

.

You can also create sine waves with orientations along any direction by using both `X`

and `Y`

as arguments in `np.sin()`

:

# exploring_meshgrid.py import matplotlib.pyplot as plt import numpy as np wavelength = 0.5 x = y = np.linspace(-1, 1, 1000) X, Y = np.meshgrid(x, y) output = np.sin(2 * np.pi * (X + 2 * Y) / wavelength) plt.imshow( output, cmap="gray", extent=[np.min(x), np.max(x), np.min(y), np.max(y)] ) plt.show()

This code gives the following sine wave:

### Another 2D Example

Let’s plot another 2D equation. You replace the sine function with a Gaussian function in the following code. You also swap to use Matplotlib’s object-oriented interface for plotting. This interface gives you more flexibility to customise your plot:

# exploring_meshgrid.py import matplotlib.pyplot as plt import numpy as np wavelength = 0.5 x = y = np.linspace(-1, 1, 1000) X, Y = np.meshgrid(x, y) output = np.exp(-5 * (X ** 2) / 2 - 5 * (Y ** 2) / 2) fig = plt.figure() ax = fig.add_subplot(121) ax.imshow( output, cmap="copper", extent=[np.min(x), np.max(x), np.min(y), np.max(y)], ) ax = fig.add_subplot(122, projection="3d") ax.plot_surface(X, Y, output, cmap="copper") plt.show()

You create two subplots within the same figure. The first one is the 2D view of the array. The second plot shows a 3D projection. The output from this code is the following:

## Using `numpy.meshgrid()`

With Higher Dimensions

As with everything else in maths and science, it gets trickier to visualise things in higher dimensions. But let’s just push one dimension further.

Let’s go back to the simpler vectors which had five values each. You can extend this further into three dimensions by adding another vector representing the third dimension. You will switch to using the `"ij"`

indexing system. For clarity, you can use *(i, j, k)* instead of *(x, y, z)*. The third vector has the values `20`

to `24`

. As I did earlier, I’m using different values for each axis for clarity:

>>> import numpy as np >>> i = np.array(range(5)) >>> i array([0, 1, 2, 3, 4]) >>> j = np.array(range(12, 17)) >>> j array([12, 13, 14, 15, 16]) >>> k = np.array(range(20, 25)) >>> k array([20, 21, 22, 23, 24])

You can use `numpy.meshgrid()`

in the same way you did for the 2D case. In this case, you need to use the optional parameter `indexing`

to use the `"ij"`

indexing system:

>>> I, J, K = np.meshgrid(i, j, k, indexing="ij") >>> I array([[[0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0]], [[1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1]], [[2, 2, 2, 2, 2], [2, 2, 2, 2, 2], [2, 2, 2, 2, 2], [2, 2, 2, 2, 2], [2, 2, 2, 2, 2]], [[3, 3, 3, 3, 3], [3, 3, 3, 3, 3], [3, 3, 3, 3, 3], [3, 3, 3, 3, 3], [3, 3, 3, 3, 3]], [[4, 4, 4, 4, 4], [4, 4, 4, 4, 4], [4, 4, 4, 4, 4], [4, 4, 4, 4, 4], [4, 4, 4, 4, 4]]]) >>> J >>> array([[[12, 12, 12, 12, 12], [13, 13, 13, 13, 13], [14, 14, 14, 14, 14], [15, 15, 15, 15, 15], [16, 16, 16, 16, 16]], [[12, 12, 12, 12, 12], [13, 13, 13, 13, 13], [14, 14, 14, 14, 14], [15, 15, 15, 15, 15], [16, 16, 16, 16, 16]], [[12, 12, 12, 12, 12], [13, 13, 13, 13, 13], [14, 14, 14, 14, 14], [15, 15, 15, 15, 15], [16, 16, 16, 16, 16]], [[12, 12, 12, 12, 12], [13, 13, 13, 13, 13], [14, 14, 14, 14, 14], [15, 15, 15, 15, 15], [16, 16, 16, 16, 16]], [[12, 12, 12, 12, 12], [13, 13, 13, 13, 13], [14, 14, 14, 14, 14], [15, 15, 15, 15, 15], [16, 16, 16, 16, 16]]]) >>> K >>> array([[[20, 21, 22, 23, 24], [20, 21, 22, 23, 24], [20, 21, 22, 23, 24], [20, 21, 22, 23, 24], [20, 21, 22, 23, 24]], [[20, 21, 22, 23, 24], [20, 21, 22, 23, 24], [20, 21, 22, 23, 24], [20, 21, 22, 23, 24], [20, 21, 22, 23, 24]], [[20, 21, 22, 23, 24], [20, 21, 22, 23, 24], [20, 21, 22, 23, 24], [20, 21, 22, 23, 24], [20, 21, 22, 23, 24]], [[20, 21, 22, 23, 24], [20, 21, 22, 23, 24], [20, 21, 22, 23, 24], [20, 21, 22, 23, 24], [20, 21, 22, 23, 24]], [[20, 21, 22, 23, 24], [20, 21, 22, 23, 24], [20, 21, 22, 23, 24], [20, 21, 22, 23, 24], [20, 21, 22, 23, 24]]])

You can compare this output with the pictorial representation below which shows the 1D vector `i`

and its 3D counterpart `I`

. I have reluctantly opted not to use my preferred hand-drawn diagrams style for this 3D representation!

This is the 1D vector `i`

:

And this is its 3D counterpart `I`

:

The output from the code above which displays `I`

shows five 2D arrays. These represent the five horizontal layers in the diagram above.

A quick note on indexing systems: the only difference between the `"xy"`

and `"ij"`

indexing systems is that the first two dimensions are swapped. In higher-dimension arrays, the rest of the dimensions are unchanged.

## Do You Need To Use `numpy.meshgrid()`

?

I mentioned earlier that I come from a MATLAB background originally and I was used to MATLAB’s meshgrid function. When I came to Python and NumPy and spotted a function with the same name, I didn’t look any further. I had found the function I needed.

However, in NumPy, you don’t need to use `meshgrid()`

. In some cases, you’re better off not using it at all.

Let’s explore some of the options. I’ll use the 2D sine wave oriented along a diagonal that you used earlier:

# exploring_meshgrid.py import matplotlib.pyplot as plt import numpy as np wavelength = 0.5 x = y = np.linspace(-1, 1, 1000) X, Y = np.meshgrid(x, y) output = np.sin(2 * np.pi * (X + 2 * Y) / wavelength) plt.imshow( output, cmap="gray", extent=[np.min(x), np.max(x), np.min(y), np.max(y)] ) plt.show()

### Using 1D Arrays

If you’re using NumPy, you probably heard the term *broadcasting*. I’ve got an article in the pipeline about this topic which will expand on this Twitter thread about broadcasting, but in the meantime, you can read the excellent NumPy docs.

You can create a 2D array from arithmetic operations with 1D arrays as long as the shapes of the arrays follow NumPy’s broadcasting rules.

Let’s look at an example with small arrays first:

>>> import numpy as np >>> x = np.array([0, 1, 2, 3]) >>> y = np.array([0, 1, 2, 3]) >>> x array([0, 1, 2, 3]) >>> y array([0, 1, 2, 3]) >>> x.shape (4,) >>> y.shape (4,)

Arithmetic operations between these two arrays will give another 1D array:

>>> x + y array([0, 2, 4, 6]) >>> x * y array([0, 1, 4, 9])

You need one of these arrays to be transposed so that its values are along the second axis. However, you’ll first need to convert these arrays into 2D arrays in which one of the dimensions is equal to `1`

:

# Notice the extra pair of [ ] # The trailing comma is not required # but can serve as a visual reminder # that this is a 2D array, helping readability >>> x = np.array([[0, 1, 2, 3],]) >>> y = np.array([[0, 1, 2, 3],]) # Again, there are two sets of [[ ]] >>> x array([[0, 1, 2, 3]]) >>> y array([[0, 1, 2, 3]]) >>> x.shape (1, 4) >>> y.shape (1, 4) # Now, you can transpose y >>> y = y.T >>> y array([[0], [1], [2], [3]]) >>> y.shape (4, 1)

Notice that you use two sets of square brackets `[[...]]`

when creating `x`

and `y`

. This indicates that the array has two dimensions instead of 1. The trailing comma is not needed. However, it can help with readability as it’s a visual reminder that there’s another dimension included.

What happens now when you perform arithmetic operations using `x`

and `y`

?

>>> x + y array([[0, 1, 2, 3], [1, 2, 3, 4], [2, 3, 4, 5], [3, 4, 5, 6]]) >>> x * y array([[0, 0, 0, 0], [0, 1, 2, 3], [0, 2, 4, 6], [0, 3, 6, 9]])

Broadcasting rules lead to a 2D array of shape `(4, 4)`

.

You can now use the same method to refactor the code creating the 2D sine wave:

# exploring_meshgrid.py import matplotlib.pyplot as plt import numpy as np wavelength = 0.5 x = np.array([np.linspace(-1, 1, 1000)],) y = x.T output = np.sin(2 * np.pi * (x + 2 * y) / wavelength) plt.imshow( output, cmap="gray", extent=[np.min(x), np.max(x), np.min(y), np.max(y)] ) plt.show()

When you create `x`

, the array returned by `linspace()`

is one of the dimensions of the final array. The other dimension is `1`

, and the shape of `x`

is `(1, 1000)`

.

You create `y`

by transposing `x`

so that `y`

‘s shape is `(1000, 1)`

. You no longer need the 2D versions `X`

and `Y`

. Therefore, now you can use the 1D arrays `x`

and `y`

in the line that creates the 2D sine wave. The output from this code is identical to the version using `meshgrid()`

.

However, you’re not using as much memory in this version since `x`

and `y`

, although technically 2D arrays, have one of their dimensions equal to 1.

`x`

and `y`

contain `1,000`

elements each whereas `X`

and `Y`

contain `1,000,000`

elements each. This means that `X`

and `Y`

also use up `1,000`

times more memory than `x`

and `y`

:

>>> x = np.array([np.linspace(-1, 1, 1000)],) >>> x.size 1000 >>> x.nbytes 8000 # X created from meshgrid in the code used earlier >>> x = y = np.linspace(-1, 1, 1000) >>> X, Y = np.meshgrid(x, y) >>> X.size 1000000 >>> X.nbytes 8000000

This will have an impact on the speed of execution, too. You can perform a quick check by using the `timeit`

module. In this version, you increase the size of the arrays further to compare the two versions using larger arrays. I’m using a new script called `exploring_meshgrid_perf.py`

for this:

# exploring_meshgrid_perf.py import numpy as np import timeit n = 10_000 wavelength = 0.5 def using_meshgrid(): x = y = np.linspace(-1, 1, n) X, Y = np.meshgrid(x, y) output = np.sin(2 * np.pi * (X + 2 * Y) / wavelength) def using_1d_arrays(): x = np.array( [ np.linspace(-1, 1, n), ] ) y = x.T output = np.sin(2 * np.pi * (x + 2 * y) / wavelength) print( f"Using meshgrid:\n" f"{timeit.timeit('using_meshgrid()', number=1, globals=globals())}" ) print( f"\nUsing 1D arrays:\n" f"{timeit.timeit('using_1d_arrays()', number=1, globals=globals())}" )

The output from this test shows that using `meshgrid()`

is roughly three times slower than using 1D arrays for this example:

Using meshgrid: 3.360328125 Using 1D arrays: 1.1643092080000002

Performance results will vary on different setups. There is another benefit of using 1D arrays instead of `meshgrid()`

. You can use larger arrays without running out of memory. On the macOS system I’m using, I get the following outputs when I change `n`

to `50_000`

in the example above. I’ll run the two versions one at a time, starting with the 1D array version:

# exploring_meshgrid_perf.py import numpy as np import timeit n = 50_000 wavelength = 0.5 def using_meshgrid(): x = y = np.linspace(-1, 1, n) X, Y = np.meshgrid(x, y) output = np.sin(2 * np.pi * (X + 2 * Y) / wavelength) def using_1d_arrays(): x = np.array( [ np.linspace(-1, 1, n), ] ) y = x.T output = np.sin(2 * np.pi * (x + 2 * y) / wavelength) # print( # f"Using meshgrid:\n" # f"{timeit.timeit('using_meshgrid()', number=1, globals=globals())}" # ) print( f"\nUsing 1D arrays:\n" f"{timeit.timeit('using_1d_arrays()', number=1, globals=globals())}" )

The 1D array version took a while to run, but I got a result in the end:

Using 1D arrays: 224.92681420899999

Next, I tried the same script for the `meshgrid()`

version:

# exploring_meshgrid_perf.py import numpy as np import timeit n = 50_000 wavelength = 0.5 def using_meshgrid(): x = y = np.linspace(-1, 1, n) X, Y = np.meshgrid(x, y) output = np.sin(2 * np.pi * (X + 2 * Y) / wavelength) def using_1d_arrays(): x = np.array( [ np.linspace(-1, 1, n), ] ) y = x.T output = np.sin(2 * np.pi * (x + 2 * y) / wavelength) print( f"Using meshgrid:\n" f"{timeit.timeit('using_meshgrid()', number=1, globals=globals())}" ) # print( # f"\nUsing 1D arrays:\n" # f"{timeit.timeit('using_1d_arrays()', number=1, globals=globals())}" # )

But this time, there was no output as I ran out of memory:

Process finished with exit code 137 (interrupted by signal 9: SIGKILL)

You can account for some of the memory issues by using the `copy`

parameter in `meshgrid()`

and setting it to `False`

. This will create a view rather than a copy. Make sure you’re comfortable with the difference between views and copies if you’d like to use this option.

But, before deciding which of these options you prefer, you can have a look at more alternatives that are available to you in the next sections.

### Using `numpy.mgrid()`

NumPy has another way of creating a mesh-grid. Instead of using the function `meshgrid()`

, which needs 1D arrays as input arguments, you can use `numpy.mgrid`

. To use `mgrid`

you can index it directly using slices as you would with any NumPy array.

Since `mgrid`

relies on indexing in the same manner as other NumPy arrays, it’s comparable to `meshgrid()`

when using the `"ij"`

indexing system:

>>> import numpy as np >>> x = np.array(range(5)) >>> y = np.array(range(12, 17)) >>> X, Y = np.meshgrid(x, y, indexing="ij") >>> X_, Y_ = np.mgrid[:5, 12:17] >>> X array([[0, 0, 0, 0, 0], [1, 1, 1, 1, 1], [2, 2, 2, 2, 2], [3, 3, 3, 3, 3], [4, 4, 4, 4, 4]]) >>> X_ array([[0, 0, 0, 0, 0], [1, 1, 1, 1, 1], [2, 2, 2, 2, 2], [3, 3, 3, 3, 3], [4, 4, 4, 4, 4]]) >>> Y array([[12, 13, 14, 15, 16], [12, 13, 14, 15, 16], [12, 13, 14, 15, 16], [12, 13, 14, 15, 16], [12, 13, 14, 15, 16]]) >>> Y_ array([[12, 13, 14, 15, 16], [12, 13, 14, 15, 16], [12, 13, 14, 15, 16], [12, 13, 14, 15, 16], [12, 13, 14, 15, 16]])

The arrays `X`

and `Y`

obtained from `meshgrid()`

are the same as `X_`

and `Y_`

from `mgrid`

.

Since `mgrid`

uses the normal indexing notation with square brackets `[ ]`

, you can use slices with the `step`

parameter, too. The following example shows a 1D `mgrid`

:

>>> grid = np.mgrid[1:20:2] >>> grid array([ 1, 3, 5, 7, 9, 11, 13, 15, 17, 19])

This is equivalent to using `np.arange(1, 20, 2)`

. What if you want to choose the number of points in the array rather than the step size? This would be equivalent to using `linspace()`

instead of `arange()`

.

`mgrid`

has the solution for that as well. If you use a complex number as the step parameter, this value is used to represent the number of points in the array:

>>> grid = np.mgrid[1:20:30j] >>> grid array([ 1. , 1.65517241, 2.31034483, 2.96551724, 3.62068966, 4.27586207, 4.93103448, 5.5862069 , 6.24137931, 6.89655172, 7.55172414, 8.20689655, 8.86206897, 9.51724138, 10.17241379, 10.82758621, 11.48275862, 12.13793103, 12.79310345, 13.44827586, 14.10344828, 14.75862069, 15.4137931 , 16.06896552, 16.72413793, 17.37931034, 18.03448276, 18.68965517, 19.34482759, 20. ])

The array now has `30`

elements since you used `30j`

as the third parameter in the slice. This is equivalent to using `np.linspace(1, 20, 30)`

You can now add this version to `exploring_meshgrid_perf.py`

to test its performance:

# exploring_meshgrid_perf.py import numpy as np import timeit n = 10_000 wavelength = 0.5 def using_meshgrid(): x = y = np.linspace(-1, 1, n) X, Y = np.meshgrid(x, y, indexing="ij") output = np.sin(2 * np.pi * (X + 2 * Y) / wavelength) def using_1d_arrays(): x = np.array( [ np.linspace(-1, 1, n), ] ) y = x.T output = np.sin(2 * np.pi * (x + 2 * y) / wavelength) def using_mgrid(): X, Y = np.mgrid[-1:1:n*1j, -1:1:n*1j] output = np.sin(2 * np.pi * (X + 2 * Y) / wavelength) print( f"Using meshgrid:\n" f"{timeit.timeit('using_meshgrid()', number=10, globals=globals())}" ) print( f"\nUsing 1D arrays:\n" f"{timeit.timeit('using_1d_arrays()', number=10, globals=globals())}" ) print( f"\nUsing mgrid:\n" f"{timeit.timeit('using_mgrid()', number=10, globals=globals())}" )

Note that the number of runs in `timeit.timeit()`

is now `10`

for all cases and `n`

is back to `10_000`

.

The results I got when I ran this on my system were the following:

Using meshgrid: 22.100569541 Using 1D arrays: 11.517313875000003 Using mgrid: 27.486098125000005

The performance when using `meshgrid()`

and `mgrid`

are similar. In this instance, `mgrid`

was a bit slower, but you should always treat results from timing code with a bit of caution. It’s not surprising that the performance is similar since the size of the arrays created by `mgrid`

is the same as those from `meshgrid()`

.

There is also no difference in memory usage between the two versions.

### Using `numpy.ogrid()`

There’s one more option we can use. `numpy.ogrid`

is similar to `numpy.mgrid`

but creates an *open* version of the mesh-grid. The easiest way to understand this is to see it in action. In an earlier section, you created `X_`

and `Y_`

using `mgrid`

:

>>> X_, Y_ = np.mgrid[:5, 12:17] >>> X_ array([[0, 0, 0, 0, 0], [1, 1, 1, 1, 1], [2, 2, 2, 2, 2], [3, 3, 3, 3, 3], [4, 4, 4, 4, 4]]) >>> Y_ array([[12, 13, 14, 15, 16], [12, 13, 14, 15, 16], [12, 13, 14, 15, 16], [12, 13, 14, 15, 16], [12, 13, 14, 15, 16]])

You can replace `mgrid`

with `ogrid`

to see the difference between the two:

>>> x_, y_ = np.ogrid[:5, 12:17] >>> x_ array([[0], [1], [2], [3], [4]]) >>> y_ array([[12, 13, 14, 15, 16]]) >>> x_.shape (5, 1) >>> y_.shape (1, 5)

`ogrid`

creates 2D arrays in which one of the dimensions is `1`

. `x_`

has shape `(5, 1)`

and `y_`

has shape `(1, 5)`

. Does this sound familiar? This is the same method you used earlier when creating 1D arrays instead of using `meshgrid()`

. NumPy broadcasting will take care of the rest.

You can verify that this method gives the same result in `exploring_meshgrid.py`

:

# exploring_meshgrid.py import matplotlib.pyplot as plt import numpy as np wavelength = 0.5 x, y = np.ogrid[-1:1:1000j, -1:1:1000j] output = np.sin(2 * np.pi * (x + 2 * y) / wavelength) plt.imshow( output, cmap="gray", extent=[np.min(x), np.max(x), np.min(y), np.max(y)] ) plt.show()

You’ll get the same 2D sine wave you got for all the other versions when you run this script.

To test its performance, you can also add `ogrid`

to `exploring_meshgrid_perf.py`

. You can probably guess how well this version performs due to its similarity with the 1D array version above:

# exploring_meshgrid_perf.py import numpy as np import timeit n = 10_000 wavelength = 0.5 def using_meshgrid(): x = y = np.linspace(-1, 1, n) X, Y = np.meshgrid(x, y, indexing="ij") output = np.sin(2 * np.pi * (X + 2 * Y) / wavelength) def using_1d_arrays(): x = np.array( [ np.linspace(-1, 1, n), ] ) y = x.T output = np.sin(2 * np.pi * (x + 2 * y) / wavelength) def using_mgrid(): X, Y = np.mgrid[-1:1:n*1j, -1:1:n*1j] output = np.sin(2 * np.pi * (X + 2 * Y) / wavelength) def using_ogrid(): x, y = np.ogrid[-1:1:n*1j, -1:1:n*1j] output = np.sin(2 * np.pi * (x + 2 * y) / wavelength) print( f"Using meshgrid:\n" f"{timeit.timeit('using_meshgrid()', number=10, globals=globals())}" ) print( f"\nUsing 1D arrays:\n" f"{timeit.timeit('using_1d_arrays()', number=10, globals=globals())}" ) print( f"\nUsing mgrid:\n" f"{timeit.timeit('using_mgrid()', number=10, globals=globals())}" ) print( f"\nUsing ogrid:\n" f"{timeit.timeit('using_ogrid()', number=10, globals=globals())}" )

The output of the script when I ran it on my system was:

Using meshgrid: 23.056696749999997 Using 1D arrays: 11.544664791000002 Using mgrid: 28.553866499999998 Using ogrid: 11.489304624999988

The 1D array and the `ogrid`

versions are identical in terms of performance in this test.

You can also replicate the output from `ogrid`

through `meshgrid()`

by using the `sparse`

parameter in `meshgrid()`

and setting it to `True`

:

>>> import numpy as np >>> x = np.array(range(5)) >>> y = np.array(range(12, 17)) >>> x_m, y_m = np.meshgrid(x, y, indexing="ij", sparse=True) >>> x_o, y_o = np.ogrid[:5, 12:17] >>> x_m array([[0], [1], [2], [3], [4]]) >>> x_o array([[0], [1], [2], [3], [4]]) >>> y_m array([[12, 13, 14, 15, 16]]) >>> y_o array([[12, 13, 14, 15, 16]])

## Final Words

I have to make a confession. I still use `numpy.meshgrid()`

most of the time. This comes from years of using MATLAB and the ‘muscle memory’ in writing algorithms using MATLAB’s meshgrid. In most instances, the memory and speed performance are not an issue for me.

However, I’d like to start moving towards using one of the more efficient options. Instinctively, I prefer the method in which you create the 1D arrays *manually*. However, using `ogrid`

does make the code more concise. I’m still not sure which one I’ll end up regularly using once I can wean myself off `numpy.meshgrid()`

.

Possibly, using the `sparse`

parameter with `meshgrid()`

is the solution that may work best for me. This method keeps the same logic I’m used to of using `meshgrid()`

, but comes with the performance improvements of using an open mesh-grid instead of a fleshed out one.

In this article, you’ve explored why we need to create mesh-grids from 1D vectors and what mesh-grids are. You went on to explore several ways of creating mesh-grids using NumPy.

Two of the options create *fleshed-out* grids:

`numpy.meshgrid()`

`numpy.mgrid`

Fleshed-out grids include all the elements needed for the grid. This means that all elements in the 2D array are included for a 2D mesh-grid.

The other options create *open* mesh-grids. In an open mesh-grid, only one dimension of the arrays is greater than `1`

. These options rely on NumPy’s broadcasting to create the N-dimensional arrays:

- 1D arrays using
`numpy.array([...],)`

`numpy.ogrid`

`numpy.meshgrid(..., sparse=True)`

The choice about which option you prefer to use is yours to make! This may depend on the specific application you’re working on.

Enjoy creating mesh-grids with `numpy.meshgrid()`

or one of its alternatives!

## Subscribe to

## The Python Coding Stack

Regular articles for the intermediate Python programmer or a beginner who wants to “read ahead”

[…] You can read more about meshgrid(), including alternatives that may be more efficient, in the article numpy.meshgrid(): How Does It Work? When Do You Need It? Are There Better Alternatives? […]