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')
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.
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, 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:
# 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.
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:
If cnRasterSmoothingOn is set True, the level (and hence the color) assigned to each cell is determined by interpolating the values of neighboring points in the data grid. If cnRasterSmoothingOn is False, ContourPlot creates a discrete raster plot: any raster cell whose center lies within the rectangular area bounded by lines halfway between each grid point (in data space) is given the color assigned to the level representing the grid point datum.
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
: