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. Still, we would like some means to visualize the data on the native grid, so skipping this for now…
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:
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:
#!/usr/bin/env python3
import plac
import ngl
import numpy
import xarray
datafile = '/global/cfs/cdirs/acme/inputdata/atm/cam/inic/homme/cami_mam3_Linoz_ne4np4_L72_c160909.nc'
gridfile = '/global/cfs/cdirs/acme/mapping/grids/ne4np4-pentagons_c100308.nc'
varname = 'Q'
# Read data
ds_data = xarray.open_dataset(datafile)
ds_grid = xarray.open_dataset(gridfile).rename({'grid_size': 'ncol'})
data = ds_data[varname].isel(time=0, lev=-1).squeeze()
lon_edges = ds_grid['grid_corner_lon']
lat_edges = ds_grid['grid_corner_lat']
# Setup the canvas
plot_format='png'
wks = ngl.open_wks(plot_format, './' + varname + '_markers')
# Make map plot
res = ngl.Resources()
res.cnFillOn = True
res.cnLinesOn = False
res.cnFillPalette = 'MPL_viridis'
res.cnFillMode = 'CellFill'
res.sfXCellBounds = lon_edges.values
res.sfYCellBounds = lat_edges.values
res.mpGeophysicalLineColor='white'
res.mpGridAndLimbOn = False
res.lbOrientation = 'horizontal'
res.lbTitleString = '%s (%s)'%(data.long_name, data.units),
res.nglFrame = False
map_plot = ngl.contour_map(wks, data.values, res)
# Add polymarkers for cell centers
polyres = ngl.Resources()
polyres.gsMarkerColor = 'black'
polyres.gsMarkerIndex = 1
poly_plot = ngl.polymarker(
wks, map_plot,
ds_grid['grid_center_lon'].values,
ds_grid['grid_center_lat'].values,
polyres
)
ngl.frame(wks)
# Close things
ngl.end()
This results in the following figure:
The black dots indicate the GLL nodes, and the color fill represents the values at each GLL node. It’s also possible to do this using Matplotlib, but we have to define our own function to draw the polygons:
#!/usr/bin/env python3
from matplotlib import pyplot
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from matplotlib.collections import PatchCollection, PolyCollection
from cartopy import crs
import numpy
import xarray
def plot_unstructured(xv, yv, data, antialiased=False, **kwargs):
"""
Plot unstructured data. xv and yv specify the x and y coordinates
of the vertices of each cell and should have shape [ni,nv] where ni
is the size of the grid (number of cells) and nv is the number of
vertices for each cell. Data contains the values of the data you
want to plot, and should have dimension [ni,]. The intent with this
function is that you will need to read in an auxillary SCRIP-format
file that describes the grid (unless this information is present in
the same file that contains the data; unlikely) as well as the file
that contains the data, and grab the cell corner information from
the SCRIP file and the data from the data file. This function will
then plot the data on the native grid by looping over each cell and
drawing patches for each. Note that this will probably be really
slow for big grids! Seems to work alright up to about ne120 or so,
but really some more clever techniques should probably be used here
(parallelism?).
NOTE: To avoid artifacts due to antialiasing, you should probably pass
antialiaseds=False to **kwargs.
"""
# Create array of cell vertices, indexed [npoints, ncorners, 2]
corners = numpy.stack([xv, yv], axis=2)
# Go back and fix corners where they wrap; we shouldn't have to do
# this with cartopy, but it seems we do...
for i in range(corners.shape[0]):
if any(corners[i,:,0] < -90) and any(corners[i,:,0] > 90):
corners[i,:,0] = numpy.where(corners[i,:,0] < -90, corners[i,:,0] + 360, corners[i,:,0])
if any(corners[i,:,1] < -45) and any(corners[i,:,1] > 45):
corners[i,:,1] = numpy.where(corners[i,:,1] < -45, corners[i,:,1] + 90, corners[i,:,1])
# Create a PatchCollection from our aggregated list of PathPatches
p = PolyCollection(corners, array=data, antialiaseds=antialiased, **kwargs)
# Add the collection to the axes
ax = pyplot.gca()
ax.add_collection(p)
# Set sane axes limits
ax.set_xlim([xv.min(), xv.max()])
ax.set_ylim([yv.min(), yv.max()])
# Return collection of patches
return p
# Function to convert longitude from 0-360 to -180 to 180
def fix_lon(lon):
return numpy.where(lon > 180, lon - 360, lon)
# Plot example
datafile = '/global/cfs/cdirs/acme/inputdata/atm/cam/inic/homme/cami_mam3_Linoz_ne4np4_L72_c160909.nc'
gridfile = '/global/cfs/cdirs/acme/mapping/grids/ne4np4-pentagons_c100308.nc'
varname = 'Q'
# Read data
ds_data = xarray.open_dataset(datafile)
ds_grid = xarray.open_dataset(gridfile).rename({'grid_size': 'ncol'})
data = ds_data[varname]
lon_edges = ds_grid['grid_corner_lon']
lat_edges = ds_grid['grid_corner_lat']
# Reduce data if we need to
if 'time' in data.dims: data = data.isel(time=0).squeeze()
if 'lev' in data.dims: data = data.isel(lev=-1).squeeze()
# Make plot
figure, ax = pyplot.subplots(
1, 1,
#subplot_kw=dict(projection=crs.Orthographic(central_longitude=-100))
subplot_kw=dict(projection=crs.PlateCarree())
)
pl = plot_unstructured(
fix_lon(lon_edges.values), lat_edges.values, data.values,
transform=crs.PlateCarree()
)
ax.set_global()
ax.coastlines(color='white', linewidth=0.5)
# Add colorbar to plot
cb = pyplot.colorbar(
pl, orientation='horizontal',
label='%s (%s)'%(data.long_name, data.units), pad=0.05
)
# Add markers to plot for GLL centers
pl = ax.plot(
fix_lon(ds_grid['grid_center_lon'].values),
ds_grid['grid_center_lat'].values,
marker=',', markersize=0.5,
MarkerFaceColor='black', MarkerEdgeColor='black',
linestyle='none',
)
# Save plot
plotfile = '%s_example.png'%varname
figure.savefig(plotfile, bbox_inches='tight')
The resulting plot:
It may be possible to simplify this example, but so far this is the best solution I’ve come up with to do this in mpl. Regardless, both of these approaches have the downside of representing the GLL point data as cell-centered data, which is maybe not completely appropriate for data on the GLL grid (however, it probably is the best representation of data on the alternative pseudo-finite-volume physics grids, e.g, ne4pg2).
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:
#!/usr/bin/env python3
from matplotlib import pyplot
from cartopy import crs
import numpy
import xarray
# Function to convert longitude from 0-360 to -180 to 180
def fix_lon(lon):
return numpy.where(lon > 180, lon - 360, lon)
datafile = '/global/cfs/cdirs/acme/inputdata/atm/cam/inic/homme/cami_mam3_Linoz_ne4np4_L72_c160909.nc'
gridfile = '/global/cfs/cdirs/acme/mapping/grids/ne4np4-pentagons_c100308.nc'
varname = 'Q'
# Read data
ds_data = xarray.open_dataset(datafile)
data = ds_data[varname]
lon = ds_data['lon']
lat = ds_data['lat']
# Reduce data if we need to
if 'time' in data.dims: data = data.isel(time=0).squeeze()
if 'lev' in data.dims: data = data.isel(lev=-1).squeeze()
# Make plot
figure, ax = pyplot.subplots(
1, 1,
#subplot_kw=dict(projection=crs.Orthographic(central_longitude=-100))
subplot_kw=dict(projection=crs.PlateCarree())
)
pl = ax.tripcolor(
fix_lon(lon.values), lat.values, data.values,
transform=crs.PlateCarree()
)
ax.set_global()
ax.coastlines(color='white', linewidth=0.5)
# Add colorbar to plot
cb = pyplot.colorbar(
pl, orientation='horizontal',
label='%s (%s)'%(data.long_name, data.units), pad=0.05
)
# Add markers to plot for GLL centers
pl = ax.plot(
fix_lon(lon.values),
lat.values,
marker=',', markersize=0.5,
MarkerFaceColor='black', MarkerEdgeColor='black',
linestyle='none',
)
# Save plot
plotfile = '%s_tri_example.png'%varname
figure.savefig(plotfile, bbox_inches='tight')
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. 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, it can be expected that most of the time is spent doing the triangulation; the plotting itself should be relatively fast (TODO: confirm). Thus, we should be able to perform the triangulation once per grid, and save the results to use at plotting time. 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:
# Read data to work with
datafile = '/global/cfs/cdirs/acme/inputdata/atm/cam/inic/homme/cami_mam3_Linoz_ne4np4_L72_c160909.nc'
scrip_file = '/global/cfs/cdirs/acme/mapping/grids/ne4np4-pentagons_c100308.nc'
latlon_file = '/global/cfs/cdirs/acme/mapping/grids/ne4np4_latlon_c100308.nc'
quad_file = '/global/cscratch1/sd/bhillma/acme_scratch/cori-knl/master.ne4_ne4.F-EAMv1-AQP1.test/run/SEMapping.nc'
datasets = {
'data': xarray.open_dataset(datafile),
'scrip': xarray.open_dataset(scrip_file),
'latlon': xarray.open_dataset(latlon_file),
'quad': xarray.open_dataset(quad_file),
}
# converts quad elements into tri elements
def quads_to_tris(quads):
tris = numpy.zeros([2 * quads.shape[0],3])
for i in range(quads.shape[0]):
j = 2 * i
tris[j,0] = quads[i,3] - 1
tris[j,1] = quads[i,2] - 1
tris[j,2] = quads[i,1] - 1
tris[j + 1,0] = quads[i,1] - 1
tris[j + 1,1] = quads[i,0] - 1
tris[j + 1,2] = quads[i,3] - 1
return tris
from matplotlib.tri import Triangulation
figure = pyplot.figure(figsize=(15, 10))
ax = pyplot.axes(projection=crs.PlateCarree())
quads = datasets['quad']['element_corners']
triangles = quads_to_tris(quads.transpose().values)
triangulation = Triangulation(fix_lon(x), y, triangles)
pl = ax.tripcolor(triangulation, data, shading='flat')
ax.set_global()
ax.coastlines(linewidth=0.2)
pl = ax.plot(fix_lon(x), y, 'k.', markersize=1)
figure.savefig('%s_garbage.png'%varname, bbox_inches='tight')
But so far, my attempts to use this to plot data that wraps around the dateline has resulted in garbage.