18

The matplotlib.pyplot.contour() function takes 3 input arrays X, Y and Z.
The arrays X and Y specify the x- and y-coordinates of points, while Z specifies the corresponding value of the function of interest evaluated at the points.

I understand that np.meshgrid() makes it easy to produce arrays which serve as arguments to contour():

X = np.arange(0,5,0.01)
Y = np.arange(0,3,0.01)

X_grid, Y_grid = np.meshgrid(X,Y)
Z_grid = X_grid**2 + Y_grid**2

plt.contour(X_grid, Y_grid, Z_grid)  # Works fine

This works fine. And conveniently, this works fine too:

plt.contour(X, Y, Z_grid)  # Works fine too

However, why is the Z input required to be a 2D-array?

Why is something like the following disallowed, even though it specifies all the same data aligned appropriately?

plt.contour(X_grid.ravel(), Y_grid.ravel(), Z_grid.ravel())  # Disallowed

Also, what are the semantics when only Z is specified (without the corresponding X and Y)?

1
  • linspace can't have a float for its third (optional) argument in your example.
    – decvalts
    Commented Feb 4, 2017 at 22:50

5 Answers 5

23

Looking at the documentation of contour one finds that there are a couple of ways to call this function, e.g. contour(Z) or contour(X,Y,Z). So you'll find that it does not require any X or Y values to be present at all.

However in order to plot a contour, the underlying grid must be known to the function. Matplotlib's contour is based on a rectangular grid. But even so, allowing contour(z), with z being a 1D array, would make it impossible to know how the field should be plotted. In the case of contour(Z) where Z is a 2D array, its shape unambiguously sets the grid for the plot.

Once that grid is known, it is rather unimportant whether optional X and Y arrays are flattened or not; which is actually what the documentation tells us:

X and Y must both be 2-D with the same shape as Z, or they must both be 1-D such that len(X) is the number of columns in Z and len(Y) is the number of rows in Z.

It is also pretty obvious that someting like plt.contour(X_grid.ravel(), Y_grid.ravel(), Z_grid.ravel()) cannot produce a contour plot, because all the information about the grid shape is lost and there is no way the contour function could know how to interprete the data. E.g. if len(Z_grid.ravel()) == 12, the underlying grid's shape could be any of (1,12), (2,6), (3,4), (4,3), (6,2), (12,1).

A possible way out could of course be to allow for 1D arrays and introduce an argument shape, like plt.contour(x,y,z, shape=(6,2)). This however is not the case, so you have to live with the fact that Z needs to be 2D.

However, if you are looking for a way to obtain a countour plot with flattened (ravelled) arrays, this is possible using plt.tricontour().

plt.tricontour(X_grid.ravel(), Y_grid.ravel(), Z_grid.ravel()) 

Here a triangular grid will be produced internally using a Delaunay Triangualation. Therefore even completely randomized points will produce a nice result, as can be seen in the following picture, where this is compared to the same random points given to contour.

enter image description here

(Here is the code to produce this picture)

5
  • Thanks for the response. I agree that the information about grid shape is lost when Z is passed as the sole argument. However, in the case of plt.contour(X_grid.ravel(), Y_grid.ravel(), Z_grid.ravel()), we simply need to examine the range of values in the first two input arrays, and that would determine the precise layout of the grid points.
    – dhrumeel
    Commented Feb 4, 2017 at 23:06
  • @dhrumeel Unfortunately, your comment is incorrect. Consider the following 12-element 1D lists: x = [-0.04 1.04 2.02 0.03 1. 1.93 0.01 1.06 1.94 -0.07 0.99 1.95]; y= [ 0.04 0.17 0.14 1.9 1.85 1.74 3.61 3.69 3.51 5.23 5.33 5.27]; z = [-0.04 0.85 0.89 -0.01 -0.23 -0.16 -0.01 -0.74 -0.87 -0.03 0.48 0.49] which shape of the grid of plt.contour(x,y,z) could you infer from that? (Those are valid points and have been produced from 2D arrays using ravel). Commented Feb 4, 2017 at 23:27
  • Ok I tried contour plots by reshaping your data to (3,4), (4,3), (6,2) etc., and I do get different results in all cases. So I see your point, but I still don't get where the 'extra' information for the contours is coming from. In all cases, we are explicitly spelling out the x and y-coordinates of the points with the corresponding height of the terrain at those points. So plt.contour() is somehow using the 'orientation' of the points in the array as well as their actual locations based on their coordinate values?
    – dhrumeel
    Commented Feb 5, 2017 at 1:02
  • Yes there is some intrinsic assumption on the order of the points on the grid. See Ilya's answer and the link to the code therein. I updated the answer to incluce tricontour which acutally accepts ravelled arrays. Commented Feb 5, 2017 at 16:40
  • @ImportanceOfBeingErnest Thank you for such an interesting answer, I have posted this question, it has been well over two weeks I am trying to get a 4D contour plot. But creating the grid for a random X,Y,Z data and then counting the data points in each grid and then getting their average seems to be a challenge. I was wondering if you could kindly take a look at my question and see if you can help me out? Many thanks!
    – Pygin
    Commented May 8, 2021 at 17:46
4

The actual code of an algorithm behind plt.contour can be found in _countour.cpp. It is rather complicated C-code, so it is difficult to follow it precisely, but if I were trying to make some contours-generating code I would do it in the following way. Pick some point (x, y) at the border and fix its z-value. Iterate over nearby points and pick that one for which z-value is the closest to z-value of the first point. Continue iteration for new point, pick nearby point with the z-value closest to the desired (but check that you do not return to a point you just visited, so you have to go in some "direction"), and continue until you get a cycle or reach some border.

It seems that something close (but a bit more complex) is implemented in _counter.cpp.

As you see from the informal description of the algorithm, to proceed you have to find a point which is "nearby" to the current one. It is easy to do if you have a rectangular grid of points (need about 4 or 8 iterations like this: (x[i+1][j], y[i+1][j]), (x[i][j+1], y[i][j+1]), (x[i-1][j], y[i-1][j]) and so on). But if you have some randomly selected points (without any particular order), this problem becomes difficult: you have to iterate over all the points you have to find nearby ones and make the next step. The complexity of such step is O(n), where n is a number of points (typically a square of a size of a picture). So an algorithm becomes much slower if you don't have a rectangular grid.

This is why you actually need three 2d-arrays that correpsponds to x's, y's and z's of some points located over some rectangular grid.

As you correctly mention, x's and y's can be 1d-arrays. In this case, the corresponding 2d-arrays are reconstructed with meshgrid. However, in this case you have to have z as 2d-array anyway.

If only z is specified, x and y are range's of appropriate lengths.

EDIT. You can try to "fake" two-dimensional x, y and z arrays in such a way that x and y does not form a rectangular grid to check if my assumptions are correct.

import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

x = np.random.uniform(-3, 3, size=10000)
y = np.random.uniform(-3, 3, size=10000)
z = x**2 + y**2
X, Y, Z = (u.reshape(100, 100) for u in (x, y, z))
plt.contour(X, Y, Z)

Incorrect result

As you see, the picture does not look like anything close to the correct graph if (x, y, z)'s are just some random points.

Now let us assume that x is sorted as a preprocessing step as @dhrummel suggests in the comments. Note that we can't sort x and y simultaniously as they are not independent (we want to preserve the same points).

x = np.random.uniform(-3, 3, size=10000)
y = np.random.uniform(-3, 3, size=10000)
z = x**2 + y**2
xyz = np.array([x, y, z]).T
x, y, z = xyz[xyz[:, 0].argsort()].T
assert (x == np.sort(x)).all()
X, Y, Z = (u.reshape(100, 100) for u in (x, y, z))
plt.contour(X, Y, Z)

x is sorted now

Again, the picture is incorrect, due to the fact that y's are not sorted (in every column) as they were if we had rectangular grid instead of some random points.

3
  • Thanks for the response. I see your point about how the grid makes it easier for contour to do the plotting. I wonder if this might be achieved by simply sorting the X and Y input arrays appropriately, as a preprocessing step.
    – dhrumeel
    Commented Feb 4, 2017 at 23:13
  • @dhrumeel, I added some examples to address this comment. In general, the problem with sorting is that you cannot sort X and Y simultaniously. However, to make algorithm work you need that i, j point of your matrix is close to i±1, j±1. To this end, you have x's and y's to be sorted in every single row/column. It is impossible to achieve this if you just have some random points, without rectangular structure. Commented Feb 5, 2017 at 9:10
  • Thanks, your examples clarify the point perfectly. So it seems like a hidden assumption of the algorithm that contour and contourf use, that the inputs are well-ordered in a grid.
    – dhrumeel
    Commented Feb 10, 2017 at 2:07
3

Imagine that you want to plot a three-dimensional graph. You have a set of x points and a set of y points. The goal is to produce a value z for each pair of x and y, or in other words you need a function f such that it generates a value of z so that z = f(x, y).

Here's a good example (taken from MathWorks):

surf

The x and y coordinates are at the bottom right and bottom left respectively. You will have a function f such that for each pair of x and y, we generate a z value. Therefore, in the code you have provide, the numpy.meshgrid call will generate two 2D arrays such that for each unique spatial location, we will observe the x and y value that are unique to that location.

For example, let's use a very small example:

In [1]: import numpy as np

In [2]: x, y = np.meshgrid(np.linspace(-1, 1, 3), np.linspace(-1, 1, 3))
In [3]: x
Out[3]:
array([[-1.,  0.,  1.],
       [-1.,  0.,  1.],
       [-1.,  0.,  1.]])

In [4]: y
Out[4]:
array([[-1., -1., -1.],
       [ 0.,  0.,  0.],
       [ 1.,  1.,  1.]])

Take a look at row number 2 and column number 1 for example (I'm starting indexing at 0 btw). This means at this spatial location, we will have coordinate x = 0. and y = 1. numpy.meshgrid gives us the x and y pair that is required to generate the value of z at that particular coordinate. It's just split up into two 2D arrays for convenience.

Now what to finally put in your z variable is that it should use the function f and process what the output is for every value in x and its corresponding y.

Explicitly, you will need to formulate a z array that is 2D such that:

z = [f(-1, -1) f(0, -1) f(1, -1)]
    [f(-1,  0) f(0,  0) f(1,  0)]
    [f(-1,  1) f(0,  1) f(1,  1)]

Look very carefully at the spatial arrangement of x and y terms. We generate 9 unique values for each pair of x and y values. The x values span from -1 to 1 and the same for y. Once you generate this 2D array for z, you can use contourf to draw out level sets so that each contour line will give you the set of all possible x and y values that equal the same value of z. In addition, in between each adjacent pair of distinct lines, we fill in the area in between by the same colour.

Let's finish this off with an actual example. Suppose we have the function f(x, y) = exp(-(x**2 + y**2) / 10). This is a 2D Gaussian with a standard deviation of sqrt(5).

Therefore, let's generate a grid of x and y values, use this to generate the z values and draw a contourf plot:

import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-1, 1, 101)
y = x
x, y = np.meshgrid(x, y)
z = np.exp(-(x**2 + y**2) / 10)       
fig,ax2 = plt.subplots(1)    
ax2.contourf(x,y,z)
plt.show()

We get:

Contour

1

The reason for X and Y to be 2D is the following. Z matches to each (x,y) coordinate in the axes system a corresponding "depth" to create a 3D plot with x,y,and z coordinates.

Now assume we want to point on a arbitrary point within the axes system. We can do that by providing the x and y coordinates (x,y) for this point.For example (0,0). Now consider the "line" with the x value 1. On this line there is a number of n y values, which looks smth like:

enter image description here

If we plot this lines for all x values and y values we will get smth. like:

enter image description here

As you can see we have a 2D annotation which is consisting of 2 2D arrays, one for the x values which has the shape:

1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 10
#--> Two dimensional x values array

and one for the y values which has the shape:

10 10 10 10 10 10 10 10 10 10
9 9 9 9 9 9 9 9 9 9 
8 8 8 8 8 8 8 8 8 8
...
1 1 1 1 1 1 1 1 1 1
0 0 0 0 0 0 0 0 0 0
#--> Two dimensional y values array

Those two together provides the (x,y) coordinates for each point within the coordinate system. Now we can plot for each point the "depth" means the Z value (z coordinate). Now it is also obvious why the Z variable must be 2 dimensional with the shape (len(x),len(y)) because otherwise it can't provide a value for all points.

This behaviour can be realised by either providing 2D x,y, and z arrays to the function OR: provide 1D x and y arrays to the function and the function internally creates the 2D mesh from the x and y values with smth. like X,Y=np.meshgrid(x,y) but nevertheless z must be two dimensional.

0

Let me explain in a simple way, since I thought Z should not be 2D as well. contourf() needs X and Y to build its own space, and the relation Z(X,Y) to build a complete space, rather than merely using several points with 1D X,Y,Z information.

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.