The purpose of this page is to document different strategies for plotting data on the spectral element cubed-sphere grid. The unstructured grid presents some challenges of it’s own, as a lot of common plotting utilities expect regular gridded data, but there are some additional concerns to plotting the spectral element data.

Regridding

Perhaps the “easiest” way to plot data on the SE grid is to not plot it on the SE grid at all, but to regrid the data first to a regular latitude-longitude grid and use standard tools to plot. This is the approach that E3SM-Diags uses, and once the data is regridded it can be quite quick to plot, outperforming the rest of the options here. It’s also pretty easy to regrid on the fly using scipy.interpolate.griddata. For example, the follow will regrid SE data pretty quickly:

from scipy.interpolate import griddata
import numpy

xi = numpy.linspace(0, 360, 2 * 360)  # to regrid to 1/2 degree
yi = numpy.linspace(-90, 90, 2 * 180)  # to regrid to 1/2 degree
data_regridded = gridddata((x, y), data, (xi[None,:], yi[:,None]), method='linear')

...

Overview

EAM is running on an unstructured finite element grid, with both spectral element data (“GLL data”: vertex data given on the GLL nodes) and “PG2 data”, finite volume physgrid” data, given as cell averages for each cell. See SE Atmosphere Grid Overview (EAM & CAM).

Regridding: Interpolate to lat/lon grid first. (both PG2 and GLL data): The most common approach is to interpolate the data to a FV latitude longitude grid, either “cap” or offset: CAM-FV Grid Overview . This allows one to use a wide variety of visualization software. Offline interpolation is quick and painless with NCO’s ncremap. Mapping files are collected on the E3SM inputdata server: https://web.lcrc.anl.gov/public/e3sm/mapping/maps/

GLL vertex data: triangulation + dual grid cell fill. (ugly): This is the standard option in NCL and pyNGL for plotting unstructured data. NCL will internally compute a Delaunay triangulation of the data (via http://www.cs.cmu.edu/~quake/triangle.html ), where the data now lives on the vertex of each triangle. NCL then appears to form the Voronoi diagram (dual grid of the triangulation). On this dual grid, the data is now located at the center of each Voronoi cells. With “RasterFill”, each of these Voronoi cells will be colored based on the value of the original data set. This results are an accurate representation of all the data points in the original field, but can be quite ugly and is a very unnatural representation of the original spectral element data.

GLL vertex data: triangulation followed + Gourard shading (best): This approach is not available in NCL, but is available in matplotlib. In this approach, a Delaunay triangulation is first computed. Then each triangle is shaded using a linear interpolating based on its three vertex values. This produces nice looking results. It also has the advantage that it is monotone and min/max/oscillation preserving. All values in the original data set will be plotted. It’s perhaps the best visual representation of GLL vertex data.

PG2 data: triangularization + Voronio cell fill (unnatural): NCL’s Delaunay approach can also be used for PG2 data: computing a Delaunay triangulation and then the associated Voronoi grid. However, this is a lot of work to construct Voronoi cells when in the PG2 case, the data is already given as a cell average over quadrilaterals. For PG2, using these quadrilaterals instead of Voronio cells (next option below) is a more accurate representation of the discrete data.

PG2 data: Cell fill using FV quads. (best for PG2 data): NCL has a CellBounds option which allows the user to specify the FV cell bounds for cell-centered data such as PG2. This option, combined with “CellFill”, each PG2 quad will be colored based on the value of the data point. This will result in a blocky plot, that correctly conveys the cell average nature of each data point. All values in the original data set will be preserved. Unfortunately this option is not directly supported with matplotlib?

Regridding

Perhaps the “easiest” way to plot data on the SE grid is to not plot it on the SE grid at all, but to regrid the data first to a regular latitude-longitude grid and use standard tools to plot. This is the approach that E3SM-Diags uses, and once the data is regridded it can be quite quick to plot, outperforming the rest of the options here. It’s also pretty easy to regrid on the fly using scipy.interpolate.griddata. For example, the follow will regrid SE data pretty quickly:

from scipy.interpolate import griddata
import numpy

xi = numpy.linspace(0, 360, 2 * 360)  # to regrid to 1/2 degree
yi = numpy.linspace(-90, 90, 2 * 180)  # to regrid to 1/2 degree
data_regridded = gridddata((x, y), data, (xi[None,:], yi[:,None]), method='linear')

This can then be plotted using the efficient pcolormesh matplotlib (mpl) function (or really any standard contouring/plotting/mapping functions), and the regridding itself is surprisingly fast. Still, we would like some means to visualize the data on the native grid, so we explore some other options below.

The algorithm above is probably using bilinear interpolation. This is acceptable, especially when upscaling. But for the best results:

For GLL data, the best algorithm to use (see Transition to TempestRemap for Atmosphere grids ) is TempestRemap’s “intbilin”, as this map is accurate, monotone, and can be used for both downscaling and upscaling. The only drawback is that it is not exactly conservative. If exact conservation is needed (so that the mass computed on the interpolation grid is the same as on the native grid), then TempestRemap’s less accurate “mono” should be used.

For PG2 data, to remap from FV (unstructured) to FV (latlon), we dont yet have an integrated bilinear option. ESMF’s “aave” option ( which can be obtained most efficiently using NCO’s ncremap) is conservative and good for downscaling, but will produce blocky plots when upscaling. For nicer looking plots when upscaling, ESMF’s “bilinear” is currently the only option.

Dual grid

Data on the SE grid is awkward to work with primary because the data on the Gauss-Lobatto-Legendre (GLL) grid exists at cell vertices (i.e., data is point values), while most existing analysis tools expect cell-centered data (i.e., data points represent values over some representative area bounded by each “cell”). We can sort of mimic this cell-centered view by constructing a new mesh that is a “finite volume dual grid” to the SE GLL grid, in which we construct bounding polygons around each GLL node and pretend that the data at the GLL node is representative of the value across the area of the polygon. When adding a new SE grid to E3SM, we construct a SCRIP-formatted file with this information, primarily because some of our remapping tools used to/still do require this cell-centered representation of the data. And example can be found in the inputdata repository for the ne4np4 grid here:

...