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.

Update 2024/7

Much of the information on this page is devoted to plotting GLL data ( given at cell vertices ) and using NCL. With the adoption of the PG2 grid in E3SM v2 (2021) and the end of further NCL development as of 2020, this GLL/NCL information is less relevant to most users.

The most relevant section in 2024 is now PG2 data: Cell fill using FV quads. (best for PG2 data), using Matplotlib’s Polycolleciton. The polycollection approach shades individual polygons (taken from the E3SM SCRIP file) based on the cell-centered data. The dealiased=False option is critical to ensure seamless boarders between polygons and fast performance.

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/

PG2 data: Cell fill using FV quads. (best for PG2 data): In this approach, each PG2 quad will be colored based on the value of the data for that cell. Matplotlib’s “polycollection” given below can do this very efficiently. NCL has a CellBounds option which allows the user to specify the FV cell bounds for cell-centered data such as PG2. At low resolutions, this approach results 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.

PG2 data: triangularization + Voronio cell fill (unnatural): NCL’s Delaunay approach (and Matplotlib’s tripcolor) 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.

GLL vertex data: triangulation followed + Gourard shading (looks good): 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. But if you zoom in, you can see the triangles, which are unnatural for GLL data. A Gourard/bilinear shading option for vertex data on quads would be a nice option and show the actual GLL grid when zoomed in.

GLL vertex data: dual grid cell fill. (unnatural for GLL data): This is the standard option in NCL and pyNGL for plotting unstructured data. NCL can 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. Alternatively, the user can read in a dual grid from a SCRIP file and specify the polygon containing each GLL node. With “RasterFill”, each of these dual grid 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 look ugly with occasional chevrons, star-like cells and triangles.

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:

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 upsampling. But for the best results:

For GLL and PG2 data, the best algorithm to use (see Recommended Mapping Procedures for E3SM Atmosphere Grids ) is TempestRemap’s “intbilin”, as this map is accurate, monotone, and can be used for both downsampling and upsampling. 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. mono is good for downsampling, but will produce blocky plots when upsampling. For nicer looking plots when upsampling, ESMF’s “bilinear” is currently the only option.

Dual grid - NCL

pyNGL can construct this dual grid from a triangulation if the user specifies sfXArray and sfYArray. We don't cover this option here, but instead cover the second option where the user reads a dual grid from the SCRIP file and specifies that via sfXCellBounds and sfYCellBounds. For PG2 data, this is the most natural approach, with the PG2 dual grid made up of uniform quadrilaterals given in the PG2 SCRIP file.

For GLL data, the dual grid approach is less natural, 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.

For the SCRIP approach, an example SCRIP file can be found in the inputdata repository for the ne4np4 grid here: https://web.lcrc.anl.gov/public/e3sm/mapping/grids/ne4np4-pentagons_c100308.nc

This grid file format contains the locations of each GLL node / cell center (grid_corner_lon, grid_corner_lat), and the vertices of the polygon constructed to represent the weight of each node. With this ancillary information, it is conceptually simple to plot the data on a map: draw the polygons for each GLL node, with color fill corresponding to the value at the GLL node consistent with a chosen colormap and scale. This appears to be what pyNGL or NCL’s contour_map function will do if the plot resources sfXCellBounds and sfYCellBounds are set to the grid_corner_lon and grid_corner_lat variables from the SCRIP-format dualgrid file. For example, we can plot the bottom-level water vapor from the ne4np4 initial condition file via the following:

This results in the following figure:

Open

The black dots indicate the GLL nodes, and the color fill represents the values at each GLL node.

Dual grid - Matplotlib

It’s also possible to do this using Matplotlib, but we have to define our own function to draw the polygons:

The resulting plot:

Open

It may be possible to simplify this example, but so far this is the best solution I’ve come up with to do this particular thing in mpl. Integration with cartopy for projections can be a little tricky, but appears to be working well as of 2024. This approach is ideal for PG2 data, but not so ideal for GLL data where the dual grid is unnatural.

Triangulation

Another option is to use a Delaunay triangulation to plot the GLL data using triangles. Matplotlib will can calculate a Delaunay triangulation, and there are some convenience plot functions to do this automatically. Plotting the data then becomes much simpler:

However, you can see that the resulting plot has represented the data by creating triangles between the GLL nodes:

Reading through the MPL documentation and source code, it appears that the default shading scheme “flat” is coloring each triangle based on the average value of the data at the three vertices that form each triangle. This is probably not what we want either, since ultimately we want to represent the GLL point data. However, there is an option to shade with a gradient fill that interpolates between each vertex. This is enabled by passing the shading=”gouraud” keyword argument to the plot function. With this change, we get:

Notice that the point at the top of South America still shows the local minimum that is evident in the cell-centered plots above, thus this method should still preserve gradients and oscillations in the data.

Note that the triangulation method can produce weird results (holes in the plot area) when the data spans a plot boundary. For example, as shown above, we need to “fix” the longitude coordinates to go from -180 to 180, or we get a big hole in the plot. This also happens with polar plots, whether we fix the longitudes first or not. A quick example to demonstrate:

results in the following two figures:

We can actually fix both of these examples by first projecting the coordinates explicitly, and then doing the triangulation, rather than doing the triangulation in lat/lon coordinates and then projecting:

The main downside of this method (triangulation) is that it becomes painfully slow for large grids. We may have a couple ways around this. First, we might expect that a significant amount of time is spent calculating the Delaunay triangulation. It turns out this is only partly true, as the following timing results show (using meshes ne4, ne30, ne120, ne240, ne512):

This shows the time to plot (or calculate a triangulation, in the case of the “triangulation” data) for the range of mesh sizes using three different methods: the dual grid approach above (mpl implementation), calling tripcolor without specifying a triangulation, and finally using tripcolor but first calculating the triangulation separately. Thus, it appears that a similar amount of time is spent in both the calculation of the triangulation and the plotting of the data (caveat: sample size of one for these timing figures, so huge uncertainty here; TODO: update with average timing over multiple passes). Still, we should be able to perform the triangulation once per grid, and save the results to use at plotting time. This should cut the actual plotting time in half, if the results above can be trusted. So the first option is to do this using the MPL routines, and just save the data. However, we already have something that looks something like this in the SEMapping.nc file that is output in each case run directory. This file contains information about the quadrilaterals that are formed by groups of four GLL nodes. Unfortunately, we don’t have a corresponding quadpcolor routine in MPL to plot data in quads directly. However, one can imagine splitting each quadrilateral in two triangles by connecting opposite vertices, and then creating a triangulation object from the resulting triangles. This is fairly easy to do in python:

But so far, my attempts to use this to plot data that wraps around the dateline has resulted in garbage.

pyNGL has a similar capability, but fill styles seem to be a little more limited (or maybe just not as well documented). Setting sfXarray and sfYArray instead of sfXCellBounds and sfYCellBounds, and choosing a fill other than CellFill will force pyNGL to do a similar triangulation. The different fill options available are AreaFill and RasterFill. RasterFill with a low resolution grid produces a similar result to the dual grid methods above, with a solid fill around each GLL node:

There is an option to “smooth” using RasterFill by setting cnRasterSmoothingOn = True. From the NCL documentation:

By default, setting this to True gives the appearance of a contour plot rather than a smoothed raster plot:

We can get back to something closer to the “gouraud” fill that mpl uses by significantly increasing the number of “contour” levels by adjusting res.cnLevelSpacingF (and also res.lbLabelStride in order to make sure labels are not on top of each other). In this case, we’ve adjusted res.cnLevelSpacingF = 0.0001 and res.lbLabelStride = 10: