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:
Code Block |
---|
#!/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:
Code Block |
---|
#!/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:
Code Block |
---|
#!/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:
Code Block |
---|
# 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.