Add function for converting MPAS meshes to viz triangles

[xylar]

Each pair of adjacent vertices on an MPAS cell is joined to the cell center to create a triangle. Interpolation weights and cell
indices are computed so that MPAS fields at cell centers can either be plotted as "flat" to show the individual cell polygons or
"gouraud" to smoothly interpolate between values from the closest cell centers.

Partially addresses #281

[xylar]

Testing

I used this function along with this script:

#!/usr/bin/env python

import xarray
import numpy
import matplotlib.pyplot as plt
from matplotlib.tri import Triangulation

from mpas_tools.viz import mesh_to_triangles


dsMesh = xarray.open_dataset('mpaso.rst.0501-01-01_00000.nc')
dsTris = mesh_to_triangles(dsMesh, periodicCopy=True)

sst = dsMesh.temperature.isel(Time=0, nVertLevels=0).values

sstTri = sst[dsTris.triCellIndices]

sstNode = (sst[dsTris.nodeCellIndices]*dsTris.nodeCellWeights).sum('nInterp')

nTriangles = dsTris.sizes['nTriangles']
tris = numpy.arange(3*nTriangles).reshape(nTriangles, 3)

lonNode = numpy.rad2deg(dsTris.lonNode.values).ravel()
latNode = numpy.rad2deg(dsTris.latNode.values).ravel()
sstNode = sstNode.values.ravel()

triangles = Triangulation(lonNode, latNode, tris)

plt.figure(1)
plt.tripcolor(triangles, sstNode, shading='gouraud')
plt.xlim([0., 360.])
plt.ylim([-90., 90.])

plt.figure(2)
plt.tripcolor(triangles, sstTri, shading='flat')
plt.xlim([0., 360.])
plt.ylim([-90., 90.])

plt.show()

The input mesh is a restart file from an E3SM run on the EC60to30 mesh, but could be any restart file. (If you have a mesh file that doesn't happen to have temperature as a variable, just substitute any field on cell centers for testing.)

The results for both an EC60to30 and a QU240wLI mesh (the latter showing the individual cells with "flat" plotting better):
EC60to30 flat:

EC60to30 gouraud:

QU240wLI flat:

QU240wLI gouraud:

[xylar]

@bradyrx, I haven't played around with using this approach (tripcolor) in cartopy yet. It seems kinda slow with the EC60to30 and may end up being prohibitively so for you at higher resolution. It seems like the plotting itself slow, not just the setup, but I haven't done timing to be sure.

I think the tripcolor and tricontourf approach will be more useful for us in MPAS-Analysis than anything using bokeh, datashader, etc. because I've had trouble getting the latter to work in command-line scripts.

[xylar]

I was able to make a little headway with datashader but I'm still a long way from seeing how it would work with cartopy and such.

#!/usr/bin/env python

import holoviews as hv
import xarray
import numpy
import pandas as pd
import datashader
import datashader.transfer_functions as tf
from datashader.utils import export_image

from mpas_tools.viz import mesh_to_triangles


hv.extension('matplotlib')
hv.output(dpi=300, fig='png')
hv.output(backend='matplotlib')

dsMesh = xarray.open_dataset('mpaso.rst.0501-01-01_00000.nc')
dsTris = mesh_to_triangles(dsMesh, periodicCopy=True)

sst = dsMesh.temperature.isel(Time=0, nVertLevels=0).values

sstTri = sst[dsTris.triCellIndices]

sstNode = (sst[dsTris.nodeCellIndices]*dsTris.nodeCellWeights).sum('nInterp')

nTriangles = dsTris.sizes['nTriangles']
tris = numpy.arange(3*nTriangles).reshape(nTriangles, 3)

lonNode = numpy.rad2deg(dsTris.lonNode.values).ravel()
latNode = numpy.rad2deg(dsTris.latNode.values).ravel()
sstNode = sstNode.values.ravel()

verts = pd.DataFrame({'x': lonNode, 'y': latNode, 'z': sstNode})
tris = pd.DataFrame({'v0': tris[:, 0], 'v1': tris[:, 1], 'v2': tris[:, 2]})
cvs = datashader.Canvas(plot_height=800, plot_width=1600,
                        x_range=(0., 360.), y_range=(-90., 90.))
img = tf.shade(cvs.trimesh(verts, tris), cmap=['lightblue', 'darkblue'],
               name='EC60to30')

export_image(img, filename='EC60to30_datashader')

[bradyrx]

That looks wonderful! What's the runtime on doing something like this as of now? I.e., is it something that can really be done interactively?

[xylar]

@bradyrx, For an EC30to60 mesh, definitely could be interactive. For higher res, it depends on how much of the lag I saw is one-time start-up cost of creating the geometry vs. per-render cost. I think it will be slow but might be feasible. Worth a try. Either way, this can also serve as the starting point for a holoviews or geoviews approach that might be more efficient if you have any luck figuring out how it works. Documentation and examples are pretty sparse....