I don't see how such an algorithm would be possible. In order to place a contour between two data points the contour level must be intermediate to them. If they are both bigger than the level the algorithm won't put a contour there. I'm going to close this but few free to request a reopen if you have a suggestion for how this would be possible, or if you like continue discussion at discourse.matplotlib.org. Thanks!

In order to place a contour between two data points the contour level must be intermediate to them. If they are both bigger than the level the algorithm won't put a contour there.

Sure, but you don't need to apply the intermediate value theorem if the intermediate value (up to FP tolerance) is already part of your data =)

In this case, the zeros are in your data set, so you don't need to infer them from a positive/negative pair:

>>> import numpy
>>> x = numpy.arange(-1, 1, 0.1)
>>> y = x.reshape(-1,1)
>>> f2 = numpy.abs(x - y)
>>> zeros = numpy.where(f2 == 0)
>>> print(list(zip(zeros[0],zeros[1])))
[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7), (8, 8), (9, 9), (10, 10), (11, 11), (12, 12), (13, 13), (14, 14), (15, 15), (16, 16), (17, 17), (18, 18), (19, 19)]

For contrast, I would not expect you to draw the contour if all of the values were strictly positive (again, up to some tolerance).

The basic algorithm is described at https://en.wikipedia.org/wiki/Marching_squares As described saddles are particularly troublesome to contour. I really don't think this is a practical limitation of the contouring algorithm for real data.

@orlitzky In your abs example, the only continuous coordinates in your domain that have z=0 are those that correspond to grid points on the diagonal, everywhere else has z>0.

z-values within quads are essentially bilinearly interpolated by the contouring algorithm. Consider a quad that lies along your diagonal of interest in the grid; two opposite corners of this quad have z=0 and the other two have z>0. Bilinear interpolation of these corner values gives z>0 at all locations in the quad except for the two z=0 corners.

You are assigning special significance to your particular diagonal of interest when it is not special at all. All quads are treated the same.

You could obtain the result you want by using tricontour as triangular grids don't suffer from the degeneracy of saddle quads. But you would have to arrange your triangular grid so that your preferred diagonal corresponds to edges of triangles.

Thanks for the extra information! For what it's worth, I've arrived here by way of two SageMath bug reports / use cases, the most compelling of which implicitly plots a complex f(x,y) = 0 using contour() with a single levels=[0]. Over the reals this usually poses no problems, but when the numbers are involved are complex, users are inspired to put a norm around everything since norm(z) = 0 if and only if z = 0, and the former involves only real numbers that the plotting infrastructure can handle. In other words, they plot abs(f(x,y)) = 0 instead of f(x,y) = 0 because the latter causes type errors.

Maybe in this case we can simply preprocess the data, though, and see if all of the z points lie above or below zero (inclusive). If so, the points where it touches zero (up to some tolerance...) should be the nodes of a polygon that describe the contour. I will play around a bit.