Structure of contour object returned from tricontourf

Community matplotlib-users
1

All,

We try to generate contour polygons from an unstructured triangular grid
stored in a netcdf file:

  import netCDF4
  import matplotlib.tri as tri

  # read data
  var = netCDF4.Dataset('filename.cdf').variables
  x = var['x'][:]
  y = var['y'][:]
  elems = var['element'][:,:]-1
  data = var['attrname'][:]

  # create triangulation
  triang = tri.Triangulation(x, y, triangles=elems)

  # generate contour polygons
  levels = numpy.linspace(0, maxlevel, num=numlevels)
  contour = plt.tricontourf(triang, data, levels=levels)

  # extract geometries
  for coll in contour.collections:
    # handle all geometries for one level
    for p in coll.get_paths():
      polys = p.to_polygons()
      ...

At this point we assume, that polys[0] is a linear ring to be interpreted as
a polygon exterior and polys[1:] are the corresponding interiors for
polys[0].

Here are our questions:

Is this assumption correct?
Is there any detailed documentation describing the structure of the returned
geometries?
Are the linear rings supposed to be correctly oriented (area > 0 for
exteriors and area < 0 for the interiors)?

Thanks!
Regards Hartmut

···

---------------
http://boost-spirit.com
http://stellar.cct.lsu.edu

2

Hello Hartmut,

In brief, the answers are no, no and yes.

In more detail, assuming polys is not empty then it will contain one or
more polygon exteriors and zero or more interiors, and they can be in any
order. Here is a simple example where polys[0] is an interior and polys[1]
an exterior:

    x = [0, 0, 1, 1, 0.5]
    y = [0, 1, 0, 1, 0.5]
    z = [0.5, 0.5, 0.5, 0.5, 0]
    triang = tri.Triangulation(x, y)
    contour = plt.tricontourf(triang, z, levels=[0.2, 0.4])

The returned geometries are purposefully not documented. They are an
'implementation detail' and not considered part of the public interface.
and as such they could change at any time and hence should not be relied
upon. Of course you can choose to access them if you wish, as I do myself
sometimes, but we make no promises about what the order of the polygons is,
or that it won't change tomorrow.

In reality the order of the polygons is chosen to be something that is easy
for both the contour generation routines to create and for the various
backends to render. If you were to look at the output generated by
contourf, you will see that it is organised differently from that produced
by tricontourf and is more like you would like it to be, i.e. one or more
groups of an exterior polygon followed by zero or more interiors. This is
historical as the contourf code dates from before all of the backends were
able to render arbitrary groups of exterior and interior polygons, and so
the contourf code has to calculate the order for the backends. When the
tricontourf code was written the backends were all able to calculate the
exterior/interior containment themselves, so there was no need for
tricontourf to do it as well.

Ian

···

On 26 October 2014 00:18, Hartmut Kaiser <[email protected]...> wrote:

At this point we assume, that polys[0] is a linear ring to be interpreted
as
a polygon exterior and polys[1:] are the corresponding interiors for
polys[0].

Here are our questions:

Is this assumption correct?
Is there any detailed documentation describing the structure of the
returned
geometries?
Are the linear rings supposed to be correctly oriented (area > 0 for
exteriors and area < 0 for the interiors)?

3

Thanks Ian! Your detailed answer is much appreciated.

As you might have already guessed, we have quite some problems creating clean geometries from the generated contour data. I have tried to put together one (reasonably) small test case illustrating at least one of our issues. I apologize for the lengthy code attached.

The two attached figures demonstrate what we see. Matplotlib.png (generated by the library) does not really look ok. Also, the attached shape.png shows that there are a lot of geometries generated which are self-intersecting (highlighted in dark blue), and we already skip polygons with less than 3 points. BTW, there are many polygons stacked with the same geometries.

Anything we can do about this?

Thanks!
Regards Hartmut

shape.png
shape.png733×568 3.25 KB

mathplotlib.png
mathplotlib.png815×615 19.8 KB

test.py (3.38 KB)

···

On 26 October 2014 00:18, Hartmut Kaiser <[email protected]...> wrote:
At this point we assume, that polys[0] is a linear ring to be interpreted
as
a polygon exterior and polys[1:] are the corresponding interiors for
polys[0].

Here are our questions:

Is this assumption correct?
Is there any detailed documentation describing the structure of the
returned
geometries?
Are the linear rings supposed to be correctly oriented (area > 0 for
exteriors and area < 0 for the interiors)?

Hello Hartmut,
In brief, the answers are no, no and yes.
In more detail, assuming polys is not empty then it will contain one or
more polygon exteriors and zero or more interiors, and they can be in any
order. Here is a simple example where polys[0] is an interior and
polys[1] an exterior:

    x = [0, 0, 1, 1, 0.5]
    y = [0, 1, 0, 1, 0.5]
    z = [0.5, 0.5, 0.5, 0.5, 0]
    triang = tri.Triangulation(x, y)
    contour = plt.tricontourf(triang, z, levels=[0.2, 0.4])
The returned geometries are purposefully not documented. They are an
'implementation detail' and not considered part of the public interface.
and as such they could change at any time and hence should not be relied
upon. Of course you can choose to access them if you wish, as I do myself
sometimes, but we make no promises about what the order of the polygons
is, or that it won't change tomorrow.
In reality the order of the polygons is chosen to be something that is
easy for both the contour generation routines to create and for the
various backends to render. If you were to look at the output generated
by contourf, you will see that it is organised differently from that
produced by tricontourf and is more like you would like it to be, i.e. one
or more groups of an exterior polygon followed by zero or more
interiors. This is historical as the contourf code dates from before all
of the backends were able to render arbitrary groups of exterior and
interior polygons, and so the contourf code has to calculate the order for
the backends. When the tricontourf code was written the backends were all
able to calculate the exterior/interior containment themselves, so there
was no need for tricontourf to do it as well.

---------------

http://stellar.cct.lsu.edu

4

Hartmut,

You are using masked arrays where you shouldn't, again. The documentation
for tricontour states that it expects z to be an array, it doesn't say
masked array. If you pass it a masked array, it will ignore the mask.
Hence you have a number of triangles that include a vertex with a z-value
of -99999; when contoured these are going to give you lots of thin polygons
that you don't want.

You need to stop using masked arrays where they are not expected. Your
triangulation should only contain triangles for which you have valid data
at all three vertices. So either remove invalid triangles from your 'ele'
array before creating the triangulation, or set a mask on the triangulation
once it has been created, e.g.

    point_mask_indices = numpy.where(z.mask)
    tri_mask = numpy.any(numpy.in1d(ele, point_mask_indices).reshape(-1,
3), axis=1)
    triang.set_mask(tri_mask)

Ian

···

On 1 November 2014 18:20, Hartmut Kaiser <[email protected]...> wrote:

Thanks Ian! Your detailed answer is much appreciated.

As you might have already guessed, we have quite some problems creating
clean geometries from the generated contour data. I have tried to put
together one (reasonably) small test case illustrating at least one of our
issues. I apologize for the lengthy code attached.

The two attached figures demonstrate what we see. Matplotlib.png
(generated by the library) does not really look ok. Also, the attached
shape.png shows that there are a lot of geometries generated which are
self-intersecting (highlighted in dark blue), and we already skip polygons
with less than 3 points. BTW, there are many polygons stacked with the same
geometries.

Anything we can do about this?

Thanks!
Regards Hartmut

5

Would it make sense to at least emit a warning when a mask is encountered. There are very few places in matplotlib where masked arrays are not allowed (I think histograms is the other spot, but I can’t remember for sure).

···

On Sun, Nov 2, 2014 at 11:10 AM, Ian Thomas <ianthomas23@…83…287…> wrote:

Matplotlib-users mailing list

[email protected]

https://lists.sourceforge.net/lists/listinfo/matplotlib-users

On 1 November 2014 18:20, Hartmut Kaiser <hartmut.kaiser@…287…> wrote:

Thanks Ian! Your detailed answer is much appreciated.

As you might have already guessed, we have quite some problems creating clean geometries from the generated contour data. I have tried to put together one (reasonably) small test case illustrating at least one of our issues. I apologize for the lengthy code attached.

The two attached figures demonstrate what we see. Matplotlib.png (generated by the library) does not really look ok. Also, the attached shape.png shows that there are a lot of geometries generated which are self-intersecting (highlighted in dark blue), and we already skip polygons with less than 3 points. BTW, there are many polygons stacked with the same geometries.

Anything we can do about this?

Thanks!

Regards Hartmut

Hartmut,

You are using masked arrays where you shouldn’t, again. The documentation for tricontour states that it expects z to be an array, it doesn’t say masked array. If you pass it a masked array, it will ignore the mask. Hence you have a number of triangles that include a vertex with a z-value of -99999; when contoured these are going to give you lots of thin polygons that you don’t want.

You need to stop using masked arrays where they are not expected. Your triangulation should only contain triangles for which you have valid data at all three vertices. So either remove invalid triangles from your ‘ele’ array before creating the triangulation, or set a mask on the triangulation once it has been created, e.g.

point_mask_indices = numpy.where(z.mask)

tri_mask = numpy.any(numpy.in1d(ele, point_mask_indices).reshape(-1, 3), axis=1)

triang.set_mask(tri_mask)

Ian

6

Ian,

You are using masked arrays where you shouldn't, again. The documentation
for tricontour states that it expects z to be an array, it doesn't say
masked array. If you pass it a masked array, it will ignore the
mask. Hence you have a number of triangles that include a vertex with a
z-value of -99999; when contoured these are going to give you lots of thin
polygons that you don't want.
You need to stop using masked arrays where they are not expected. Your
triangulation should only contain triangles for which you have valid data
at all three vertices. So either remove invalid triangles from your 'ele'
array before creating the triangulation, or set a mask on the
triangulation once it has been created, e.g.

    point_mask_indices = numpy.where(z.mask)
    tri_mask = numpy.any(numpy.in1d(ele, point_mask_indices).reshape(-1,
3), axis=1)
    triang.set_mask(tri_mask)

Thanks very much for this explanation! With your code everything fell into place and we see the geometries we expect to see. Great library!

Regards Hartmut

···

---------------

http://stellar.cct.lsu.edu

7

Ben,

That's certainly an option, but I'm happy with leaving it as it is with the
onus on users to read the documentation.

Ian

···

On 2 November 2014 16:30, Benjamin Root <[email protected]...> wrote:

Would it make sense to at least emit a warning when a mask is encountered.
There are very few places in matplotlib where masked arrays are not allowed
(I think histograms is the other spot, but I can't remember for sure).