mpl.tri.Triangulation does not work but plt.tricontourf does

[keltonhalbert]

As described in the title, I can successfully plot data using tricontourf but cannot successfully do it by first creating triangles using tri.Triangulation.

The reason I would want to do this is because my data is on an unstructured grid that remains static, so it would be much more efficient if I could save the triangles to a file and load them any time I need to make a plot.

Here's the error I get from Triangulation using the verbose flag:

QH6154 qhull precision error: initial facet 1 is coplanar with the interior point
ERRONEOUS FACET:
- f1
    - flags: bottom simplicial flipped
    - normal:    0.3823   0.3824  -0.8412
    - offset:   -2443442
    - vertices: p3991443(v2) p946265(v1) p3857645(v0)
    - neighboring facets: f2 f3 f4

While executing:  | qhull d Qt Qbb Qc Qz
Options selected for Qhull 2012.1 2012/02/18:
  run-id 777401308  delaunay  Qtriangulate  Qbbound-last  Qcoplanar-keep
  Qz-infinity-point  _pre-merge  _zero-centrum  Qinterior-keep  Pgood
  _max-width 1e+30  Error-roundoff 1.4e+15  _one-merge 9.7e+15
  Visible-distance 2.8e+15  U-coplanar-distance 2.8e+15  Width-outside 5.5e+15
  _wide-facet 1.7e+16

precision problems (corrected unless 'Q0' or an error)
      2 flipped facets

The input to qhull appears to be less than 3 dimensional, or a
computation has overflowed.

Qhull could not construct a clearly convex simplex from points:
- p3867215(v3): 1.9e+06 1.9e+06 -7.7e-34
- p3991443(v2): 6.4e+06    93 1.5e-17
- p946265(v1): 1e+30 1e+30 9.1e+29
- p3857645(v0):    89 6.4e+06 1.5e-17

The center point is coplanar with a facet, or a vertex is coplanar
with a neighboring facet.  The maximum round off error for
computing distances is 1.4e+15.  The center point, facets and distances
to the center point are as follows:

center point  2.5e+29  2.5e+29 2.273e+29

facet p3991443 p946265 p3857645 distance=    0
facet p3867215 p946265 p3857645 distance=    0
facet p3867215 p3991443 p3857645 distance= -2.3e+29
facet p3867215 p3991443 p946265 distance= -3.5e+13

These points either have a maximum or minimum x-coordinate, or
they maximize the determinant for k coordinates.  Trial points
are first selected from points that maximize a coordinate.

The min and max coordinates for each dimension are:
  0:     89.26     1e+30  difference= 1e+30
  1:     92.64     1e+30  difference= 1e+30
  2:  -7.704e-34     1e+30  difference= 1e+30

If the input should be full dimensional, you have several options that
may determine an initial simplex:
  - use 'QJ'  to joggle the input and make it full dimensional
  - use 'QbB' to scale the points to the unit cube
  - use 'QR0' to randomly rotate the input for different maximum points
  - use 'Qs'  to search all points for the initial simplex
  - use 'En'  to specify a maximum roundoff error less than 1.4e+15.
  - trace execution with 'T3' to see the determinant for each point.

If the input is lower dimensional:
  - use 'QJ' to joggle the input and make it full dimensional
  - use 'Qbk:0Bk:0' to delete coordinate k from the input.  You should
    pick the coordinate with the least range.  The hull will have the
    correct topology.
  - determine the flat containing the points, rotate the points
    into a coordinate plane, and delete the other coordinates.
  - add one or more points to make the input full dimensional.
Traceback (most recent call last):
  File "mpas.py", line 23, in <module>
    triang = tri.Triangulation(x, y, mask=None)
  File "/glade/apps/opt/matplotlib/1.4.3/gnu-westmere/4.8.2/lib/python2.7/site-packages/matplotlib-1.4.3-py2.7-linux-x86_64.egg/matplotlib/tri/triangulation.py", line 55, in __init__
    self.triangles, self._neighbors = _qhull.delaunay(x, y)
RuntimeError: Error in qhull Delaunay triangulation calculation: singular input data (exitcode=2)

[ianthomas23]

According to the error message, your x-coordinates span the range 89.26 to 1e30, and your y-coordinates span the range 92.64 to 1e30. Hence your domain spans 28 orders of magnitude in both dimensions. This is likely to lead to precision errors!

Why are you specifying a point at (1e30, 1e30)?

[WeatherGod]

I think the question is more of what is tricontourf() doing that
hides/fixes this problem that Triangulation isn't doing?

On Wed, May 13, 2015 at 3:08 AM, Ian Thomas notifications@github.com
wrote:

According to the error message, your x-coordinates span the range 89.26 to
1e30, and your y-coordinates span the range 92.64 to 1e30. Hence your
domain spans 28 orders of magnitude in both dimensions. This is likely to
lead to precision errors!

Why are you specifying a point at (1e30, 1e30)?

—
Reply to this email directly or view it on GitHub
#4424 (comment)
.

[keltonhalbert]

To answer the question about the 1e30 values...

I'm taking data that represents a mesh with centers having values on a lat/lon grid. I first pass through my points to Baseman to convert them to projection space, and then construct the mesh on the projection. My understanding is that 1e30 is the default masked value for data outside of the projection?

[keltonhalbert]

After looking around the tricontour and tricontourf code, in addition to the Basemap code that is working for generating the triangles, it looks like it masks the arrays before passing it to Triangulation. Perhaps this could be part of the problem?

[keltonhalbert]

After further investigation, it appears that this is in fact not a bug with Triangulation, but my poor understanding of how to handle the Basemap projection conversion.

I just copied and modified this block of code from the Basemap contourf function:

x,y = m(lons, lats)
mask = np.logical_or(x<1.e20,y<1.e20)
x = np.compress(mask,x)
y = np.compress(mask,y)