Problem about plot_trisurf of matplotlib 2.0.2

[zxgdll]

After updating matplotlib version from 1.3.1 to 2.0.2, when I want to use plot_trisurf to generate a TIN by 3d-points, I got an incomprehensible result.
I just want to know the difference of plot_trisurf between matplotlib 2.0.2 and matplotlib 1.3.1. And how i can get the similar result?
the details about this issue are referred to https://stackoverflow.com/questions/44244297/plot-trisurf-of-matplotlib-2-0-2

[ImportanceOfBeingErnest]

To make this issue a bit clearer: OP has some x,y,z data in 3 columns and wants to plot a trisurf plot. The original file has ~38000 entries. It can be found here.

To make the issue apparent, we can use a 2D tripcolor plot as well. The following script plots a scatter and a tripcolor of the points. While the scatter shows all points, the tripcolor seems to take a subset of them (or interpolate). Of course we would expect the tripcolor to look similar to the scatter and show all the triangles between points.

import matplotlib.pyplot as plt
import numpy as np
plt.rcParams["axes.xmargin"] = 0
plt.rcParams["axes.ymargin"] = 0

x,y,z = np.loadtxt('test.xyz', unpack=True)  #test.xyz contains 38010 points

fig, (ax, ax2) = plt.subplots(ncols=2)
ax.set_title("scatter")
ax2.set_title("tripcolor")

s1 = ax.scatter(x,y, c=z, s=1)
s2 = ax2.tripcolor(x,y,z)
print len(s1.get_paths()), len(s1.get_array()) #prints 1 38010
print len(s2.get_paths()), len(s2.get_array()) #prints 326 326
fig.tight_layout()
plt.show()

There seems to be no argument one can pass to the Triangulation which would help setting the resolution (like e.g. cstride for a plot_surface plot or so).

[WeatherGod]

attn @ianthomas23 , I dug through the git logs, but I am not seeing anything all that relevant. Thoughts?

[QuLogic]

I cannot get this to run with the matplotlib 1.3.1 package in conda, or one built from source:

~/code/conda/envs/mpl-old/lib/python2.7/site-packages/matplotlib/delaunay/triangulate.py:104: DuplicatePointWarning: Input data contains duplicate x,y points; some values are ignored.
  DuplicatePointWarning,
Traceback (most recent call last):
  File "issue8682.py", line 13, in <module>
    s2 = ax2.tripcolor(x,y,z)
  File "~/code/conda/envs/mpl-old/lib/python2.7/site-packages/matplotlib/axes.py", line 9195, in tripcolor
    return mtri.tripcolor(self, *args, **kwargs)
  File "~/code/conda/envs/mpl-old/lib/python2.7/site-packages/matplotlib/tri/tripcolor.py", line 62, in tripcolor
    tri, args, kwargs = Triangulation.get_from_args_and_kwargs(*args, **kwargs)
  File "~/code/conda/envs/mpl-old/lib/python2.7/site-packages/matplotlib/tri/triangulation.py", line 168, in get_from_args_and_kwargs
    triangulation = Triangulation(x, y, triangles, mask)
  File "~/code/conda/envs/mpl-old/lib/python2.7/site-packages/matplotlib/tri/triangulation.py", line 55, in __init__
    dt = delaunay.Triangulation(self.x, self.y)
  File "~/code/conda/envs/mpl-old/lib/python2.7/site-packages/matplotlib/delaunay/triangulate.py", line 124, in __init__
    self.hull = self._compute_convex_hull()
  File "~/code/conda/envs/mpl-old/lib/python2.7/site-packages/matplotlib/delaunay/triangulate.py", line 160, in _compute_convex_hull
    hull.append(edges.pop(hull[-1]))
KeyError: 31614

[zxgdll]

@QuLogic Yes, if you run with the matplotlib 1.3.1 package with this test code, you will get above error. I think it was a bug in matplotlib 1.3.1, so it was updated. If you want to see the result with matplotlib 1.3.1, you need to change s2 = ax2.tripcolor(x,y,z) to s2 = ax2.tripcolor(x[:10000],y[:10000],z[:10000]), and I hope that you can upload the result image.

[QuLogic]

Bisect puts this change to be a very long while ago:

a65066d72411e682043652afd1e51d404b85ee68 is the first bad commit
commit a65066d72411e682043652afd1e51d404b85ee68
Author: Ian Thomas <ianthomas23@gmail.com>
Date:   Mon Oct 7 21:42:28 2013 +0100

    Use qhull for delaunay triangulation

:100644 100644 509eb2aefe215f86e4ac034e992872c9f150c048 7da2e125841fbca280b15680fc9d979d237ff3ec M	CHANGELOG
:100644 100644 542205f1f0a3ae750bc08d55d9264117abaea808 92f11d0dd7fe092eccca9f46270456909bc4e0a5 M	MANIFEST.in
:040000 040000 3e626afec08d2dabe956c106ad81db2245cc89c9 33ffaf5f04401279231b516504a8c4f89255d5d3 M	examples
:040000 040000 dac54b08e0f279c7ca2193bcf169b2a024f30bc8 6fc9a566195af65d25146da387f95dac390bcfda M	lib
:100644 100644 32f4c341756cd1505ac159648222b9232a27eb8e a14e10a6da927e6f67ec5e79f02c900acf7fe510 M	setup.py
:100644 100644 c71c4d7dee2926b222ca3e874a9e017cbba26fb8 35622390e9e34942b1b06b1c90d7187e57a083fe M	setupext.py
:040000 040000 975bedd85139a8581c2c7a242ec6401642edcfa4 0353ccc0a6a00561b8efda6569b66384abdd4722 M	src

which is #2504.

[ianthomas23]

The problem is one of finite precision when calculating the Delaunay triangulation in qhull, which considers points that are near (according to a complicated definition of the word 'near') to be identical, so the triangulation is simpler than is desired. The dataset is an extreme one (in a bad way) for finite precision in that the spread of points about their mean is small (x.mean()=512767, x.max()-x.min()=134, y.mean()=303, y.max()-y.min()=5403707).

The workaround is to subtract the mean of the x, y coordinates before triangulating:

import matplotlib.tri as mtri
import numpy as np
x,y,z = np.loadtxt('test.xyz', unpack=True)
print('Number of triangles:', len(mtri.Triangulation(x, y).triangles))
print('Number of triangles after subtracting mean:', len(mtri.Triangulation(x-x.mean(), y-y.mean()).triangles))

which gives

Number of triangles: 326
Number of triangles after subtracting mean: 75840

There are lots of options in qhull for dealing with things in a different way. We do not expose those options as we don't expect our users to be qhull experts. scipy.spatial also uses qhull and does expose those options, e.g.

import scipy.spatial
import numpy as np
x,y,z = np.loadtxt('test.xyz', unpack=True)
xy = np.stack((x, y), axis=1)
triang = scipy.spatial.Delaunay(xy)
print('Number of triangles', len(triang.simplices))
triang = scipy.spatial.Delaunay(xy, qhull_options='QbB')
print('Number of triangles', len(triang.simplices))

We can conclude that the qhull options that we use are suboptimal in this case. We can change the options that we use but this is non-trivial as we need to ensure that we're not making lots of other triangulations worse to make this one better. I won't be rushing into this.

[ImportanceOfBeingErnest]

So to sum this up, this not actually an issue between different matplotlib versions and the current version is good enough to cope with the majority of use cases.
The example from above can be easily corrected by subtracting the mean from the coordinates,

s2 = ax2.tripcolor(x-x.mean(),y-y.mean(),z)

producing the desired image

The axis ticklabes can be easily corrected for. Alternatively, one may calculate the triangluation outside of matplotlib,

import scipy.spatial
xy = np.stack((x, y), axis=1)
triang = scipy.spatial.Delaunay(xy, qhull_options='QbB')
s2 = ax2.tripcolor(x,y,triang.simplices ,z)

which would yield a very similar plot as above.

I think this is good enough as an explanation and I wouldn't consider there to be any need for changes in the code.

It may be worthwhile thinking about adding a sentence like

"In some rare cases the automatic parameters for the triangulation chosen by matplotlib may not be optimal (e.g. if the mean of the points is very large compared to the distance between points). In such cases adapting the input data to center around zero or calculating the triangles with different options (see scipy.spatial.Delaunay) may be necessary."

to the documentation of tripcolor and plot_trisurface. The question is, if the number of people getting confused by this sentence wouldn't be larger than the number of people that would run into such an issue. ;-)

[zxgdll]

Thanks a lot for all replies!
This is good enough as an explanation.I will close this issue.

[ianthomas23]

@ImportanceOfBeingErnest: Your summary is only partly correct.

Old versions of mpl do give a different result. Before qhull, we used our own code for calculating Delaunay triangulations and it was pretty poor code, failing in many situations including regularly spaced grids. We will never return to using that code, qhull is here to stay.

You should not advocate using

s2 = ax2.tripcolor(x-x.mean(),y-y.mean(),z)

and then correcting the tick labels. You should remove the means to calculate the triangulation, but use the original values for rendering, e.g.

triang = mtri.Triangulation(x - x.mean(), y - y.mean())
s2 = ax2.tripcolor(x, y, triang.triangles, z)

You should not advocate people using scipy.spatial instead of matplotlib.tri so that they can use qhull options unless they understand the qhull options. The option I gave in my previous post is just an example showing qhull producing different output. It may not be the best option in this or other cases. Advocating using that option is pre-empting the work I have to do to determine the best options for mpl to cover all of our typical use cases. I am very happy for someone else to do that work, but don't tell people that it is a good option until you know it to be a good option.

We very much need to change our code, but as I said before this is non-trivial and I won't be rushing it. I won't be adding any sentence saying that the triangulation in mpl is suboptimal, I will fix the problem.

Obviously I have reopened the issue.

[tacaswell]

@ianthomas23 I think this is closed by #8873 ?

[ianthomas23]

@tacaswell Yes, it is closed by #8873.