Problem

Relevance:

The different flavors of the ICON model use a icosahedral-triangular Arakawa C grid which has most variables located in the center of each triangular cell. See here for an example visualization of a coarse global ICON grid.
As of my knowledge most of the people in the atmospheric modeling community rely on interpolation of the unstructured triangle-centered grid to a regular rectilinear lat-lon grid and do not plot the original triangular data.

Problem Description:

The current functions tricontour/tricontourf require Z to be located on the triangle vertices (x,y) and don't allow triangle-centered data:
tricontour(x, y, Z, levels, triangles=triangles, mask=mask)
tricontourf(x, y, Z, levels, triangles=triangles, mask=mask)
with x,y: (nvertices,) Z: (nvertices,) triangles: (ntri,3) mask: (ntri,)

Now analog to the implemented tripcolor functionality (see #811) allowing triangle-centered data means one of the following to be accepted:
x,y: (nvertices,) Z: (ntri,) triangles: (ntri,3) mask: (ntri,)
or alternatively taking the triangle center locations and optionally maybe neighbor triangle indices for faster computation:
x,y: (ntri,) Z: (ntri,) mask: (ntri,) [neighbors: (ntri,3,)]

One currently possible workaround (but not for masked data) is by letting matplotlib compute the delaunay triangulation of the polygon dual grid which means taking the triangle-centered data coordinates as new vertices of the dual grid. Therefore the following usage is a current workaround:
x,y: (ntri,) Z: (ntri,) triangles: (ntridelaunay,3) mask: (ntridelaunay,)
where for large grids of regular triangles: ntridelaunay ≈ 2 * ntri

Masking triangle-centered data is not possible this way because the mask from the data comes in shape (ntri,) and not in the accepted (ntridelaunay,) and alternatively np.nan values are not allowed in Z.

Example plot of tricontourf used with this workaround way with fake data values located on the centers of the yellow triangular grid (=vertices of the black hexagonal dual grid which the delaunay triangulation divides into 4 triangles per hexagon):

Proposed solution

I understood from here how the linear interpolation algorithm along the edges (black in the figure) works but as I don't understand how the currently implemented masking works I can't propose a solution to this, @ianthomas23 maybe you can explain it to us and say something about possibilities. The idea is, not to just exclude the masked data points before delaunay triangulation but to have the contour lines end at the boundaries to the masked triangles (yellow in the figure).