A simple plotting example for unstructured grid

Output from MUSICAv0 is on an unstructured grid (an unsorted vector of columns for defined latitude-longitude coordinates) making it difficult to visualize in standard plotting frameworks. NCAR has developed a small set of routines in Python as Jupyter notebooks that can be used for plotting results. A simple plotting routine with some examples is shown here:

#This first section of code is defining the modules we will be using and how they will be referenced in the code
import xarray as xr
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature

#This section of code is defining the location of the result file(s) and extracting the relevant data
result_dir = "/glade/p/acom/MUSICA/results/opt_se_cslam_nov2019.FCHIST.ne0CONUSne30x8_ne0CONUSne30x8_mt12_pecount2700_SEACRS_TS1_CAMS_final_2013_branch/atm/hist/"
case_dir = "opt_se_cslam_nov2019.FCHIST.ne0CONUSne30x8_ne0CONUSne30x8_mt12_pecount2700_SEACRS_TS1_CAMS_final_2013_branch/atm/hist/"
result_file = "opt_se_cslam_nov2019.FCHIST.ne0CONUSne30x8_ne0CONUSne30x8_mt12_pecount2700_SEACRS_TS1_CAMS_final_2013_branch.cam.h0.2013-07.nc"
#here we are using xarray to read the defined netcdf file and writing into the variable nc
nc = xr.open_dataset(result_dir+case_dir+result_file)
#we can operate on nc now, so we are taking a slice at time 0 (first entry) and 32nd level (surface)
nc1 = nc.isel(time=0,lev=31)
#from this slice we will read the latitude (lat), longitude (lon), and ozone mixing ratio (O3)
lat = nc1['lat']
lon = nc1['lon']
var_mus = nc1['O3']

#To put data on a traditional fixed-grid we are using griddata
from scipy.interpolate import griddata
#this will define the lat/lon range we are using to correspond to 0.1x0.1 over CONUS
x = np.linspace(225,300,751)
y =  np.linspace(10,60,501)
#this will put lat and lon into arrays, something that is needed for plotting
X, Y = np.meshgrid(x,y)
#here we are putting the unstructured data onto our defined 0.1x0.1 grid using linear interpolation
grid_var = griddata((lon,lat), var_mus, (X, Y), method='linear')

#We can also put the previous statements into a function that can be called later in the notebook. This simplifies the variable selection and regridding.
def grid_musica(var, time, lev):
    nctmp = nc.isel(time=time, lev=lev)
    vartmp = nctmp[var]
    X, Y = np.meshgrid(x,y)
    grid_var = griddata((lon,lat), vartmp, (X, Y), method='linear')
    return grid_var

#To simplify plotting routines we are going to define a plotting function here
def plot_map_LCC(data, lat, lon, llim, ulim, units):
    #if you want to change the figure size of map projection you can do it here
    fig, ax = plt.subplots(figsize=(16, 9), subplot_kw=dict(projection=ccrs.LambertConformal()))
    #this is the actual plotting routine, using pcolormesh. Most variables are inputs from the function call
    im = ax.pcolormesh(lon,lat,data,cmap='Spectral_r',vmin=llim, vmax=ulim, transform=ccrs.PlateCarree())
    #this is for visualization of administrative boundaries. Many shapes included in the 'Natural Earth' library can also be used
    ax.add_feature(cfeature.STATES, edgecolor='#D3D3D3')
    cb = fig.colorbar(im, ax=ax)
	return fig, ax

#Here you can see how the plotting function above is called in the notebook 
fig, ax = plot_map_LCC(grid_var,Y,X,1e-8,1e-7,'mol/mol')
#we will zoom in slightly on the extents that are defined based on map projection
ax.set_extent([-32e5, 32e5, -25e5, 22e5], crs=ccrs.LambertConformal())
plt.title("Surface O$_{3}$ July 2013 MUSICAv0 Uniform Grid")

#this section shows how easy it can be to plot new variables/domains with functions 
grid_var = grid_musica('CO',0,20)
fig, ax = plot_map_LCC(grid_var,Y,X,5e-8,1e-7,'mol/mol')
ax.set_extent([-32e5, 32e5, -25e5, 22e5], crs=ccrs.LambertConformal())
plt.title("525mb CO July 2013 MUSICAv0 Uniform Grid")
grid_var = grid_musica('CO',0,31)
fig, ax = plot_map_LCC(grid_var,Y,X,5e-8,2.5e-7,'mol/mol')
ax.set_extent([-32e5, 32e5, -25e5, 22e5], crs=ccrs.LambertConformal())
plt.title("Surface CO July 2013 MUSICAv0 Uniform Grid")

Using a pre-defined python module for plotting

For people who are not familiar with python, we have developed a python module that people can easily import and plot CAM-chem results (either structured or unstructured) with a single line command. First, download this python module for plotting (you may need to "right click" on the link and choose "Save link as..." from the menu, or find it here: https://github.com/NCAR/CAM-chem/blob/main/docs_sphinx/examples/functions/Plot_2D.py) and put it in your directory. If you are new to python, you need to know how to import a library, and you have to install required libraries for plotting (xarray, matplotlib, cartopy, etc.). See this python tutorial (or contents on this webpage) if you are new to python. 

If you want reproduce the example plots below, download these files: 

  1. structured grid (0.95 x 1.25) model output (sample.nc)
  2. unstructured grid (ne30x16) model output (sample_se.nc)
  3. SCRIP grid file for unstructured grid (sample_se_scrip.nc)

Of course, you can use your own CAM-chem history files instead of example files, unless you want to get the exactly same example plots below for checking. 

NOTE: As of 10-MAR-2021, the Plot_2D script will be maintained in CAM-chem python Github. We shall keep examples here for reference, as those could still be helpful for beginners. 

NOTE 2: As of 04-OCT-2021, the Plot_2D script can be installed using the PyPI package: https://pypi.org/project/vivaldi-a/


First, you need to import built-in python libraries and Plot_2D module you just downloaded. Make sure Plot_2D.py file is in your directory. 

Import python libraries
import xarray as xr # To read NetCDF files
from Plot_2D import Plot_2D # To draw plots
import matplotlib.cm as cm # To change colormap used in plots

Let's read the model output files using the xarray module.

Read model output files
# FV file (structured grid)
fv_filename = "sample.nc"
ds_fv = xr.open_dataset( fv_filename )

# SE with regional refinement file (unstructured grid)
se_filename = "sample_se.nc"
ds_se = xr.open_dataset( se_filename )

# scrip filename for SE with regional refinement (contains grid information)
scrip_filename = "sample_se_scrip.nc"

As a side note, you can check the variables in the file like below. 

Plotting examples

To demonstrate the examples below, you will need to use Jupyter notebook. If you are using the local machine, Jupyter notebook is OK to use (https://jupyter.org/).  However, if you use a separate server that needs remote access, you may want to install JupyterLab (https://jupyterlab.readthedocs.io/en/stable/) to use Jupyter notebook. You can still use Jupyter notebook on your server but it could be slow. If you are using Cheyenne or Casper machines, Jupyterhub is another option you can use (https://www2.cisl.ucar.edu/resources/jupyterhub-ncar). Jupyterhub is the easiest option you can start with if you are new to python, because you don't have to install python and required packages (already installed). 

A simple plot

Just pass your surface CO. "-1" corresponds to the surface level as CESM writes the output from the top of the atmosphere to the surface.

Changing colormap

You can simply change the colormap with the "cmap" keyword. Check out this website for built-in colormaps in matplotlib python library.

Multiplying scale factors

You can simply multiply any numbers (1e9 in the example below) in the input value. 

Pretty tick option

If you set pretty_tick=True, the routine will automatically find the colorbar ticks like below. 

Plotting a region of interest only

You can set "lon_range" and "lat_range" keywords to plot a region. 

Setting maximum value for the plot

Use "cmax" keyword to force the maximum value to what you want. 

Adding state/province boundary lines

You can also easily add state boundary lines with the "state" keyword.

Adding a title in the plot

You can also add a title in the plot with a title keyword. 

Plotting on unstructured grid

You can use exactly the same function and interface to plot unstructured grid model output. 

The only difference is you have to specify the location of corresponding SCRIP grid file. 

Adding a unit

You can add a unit with a unit keyword. 

Adjusting unit position can also be done using the "unit_offset" keyword.

For advanced users: Adding more custom things you want

If you are an experienced user, you can add anything you want by getting an instance. As an example, we provide how to add a location point on top of the plot. Note "p = Plot_2D(...)" instead of "Plot_2D(...)", it creates a new instance of the class and assigns this object to the local variable "p".

More custom options and keywords

There are a lot of custom options you can explore, such as changing map projections, grid line settings, coast/country/state lines, passing parent axes for multi-plot, colorbar orientation and size, etc. For more information, type "Plot_2D?" in your python or jupyter notebook shell prompt. You will find detailed information of each keyword like below. 

Other resources and tools

  • No labels


  1. I'm getting a server error trying to download the Plot_2D.py python module. Can someone look into this? Thanks.

  2. Try:  "right click" on the link and choose "Save link as..." from the menu

    Or find it here: https://github.com/NCAR/CAM-chem/blob/main/docs_sphinx/examples/functions/Plot_2D.py

  3. The right click to save file as failed for me as well. But, getting it from the github repository is even better. Thanks for that!

  4. I’ve recently begun to run MUSICA simulations on Cheyenne. I want to plot the output with Python, but have very little experience with it. I had been following this tutorial (https://github.com/NCAR/ACOM-Python-Tutorial) to try to get the JupyterLab set up, but am having issues with the environment. I contacted the developer of this tutorial and he says that it is likely out of date because so much has changed. Are there any other tutorials or resources you would recommend for getting started so that I may process the MUSICA outputs with Python? On the “Plot Output with Python” tutorial it says we have to install the required libraries for plotting (xarray, matplotlib, cartopy, etc.), so I just want to make sure I do this correctly. Thanks!

  5. CISL has been presenting a tutorial which you might find useful:


    And see the current information about JupyterHub on cheyenne/casper:


  6. The easiest way to use Jupyter notebook (and Plotting) is using JupyterHub, as Louisa mentioned above. I will include it in the main text. Thanks!

  7. Thank you both so much. I will check out these resources. 

  8. Hello,

    When I attempt create the simple plot using – Plot_2D( ds_fv['CO'][0,-1,:,:] ) – I get an index error. In my directory, I have the Plot_2D.py and sample.nc files, and my conda environment is set up. Hope you can assist me with this issue. Thanks!

    This here is my Jupyter Notebook output when i run Plot_2D( ds_fv['CO'][0,-1,:,:] ) : 

    IndexError                                Traceback (most recent call last)
    <ipython-input-5-61b9502f27a1> in <module>
    ----> 1 Plot_2D( ds_fv['CO'][0,-1,:,:] )
    ~/plot_python/sample_unstructured/Plot_2D.py in __init__(self, var, lons, lats, lon_range, lat_range, scrip_file, ax, cmap, projection, grid_line, grid_line_lw, coast, country, state, resolution, feature_line_lw, feature_color, lonlat_line, lon_interval, lat_interval, font_family, label_size, colorbar, log_scale, log_scale_min, diff, orientation, shrink, pad, fraction, extend, colorticks, colorlabels, pretty_tick, nticks, cmax, cmin, title, title_size, title_bold, unit, unit_size, unit_bold, unit_italic, unit_offset, verbose)
        555                     else:
        556                         self.cbprop = get_cbar_prop( [self.var_slice], Ntick_set=self.nticks,
    --> 557                                                     **self.kwd_pretty_tick )
        558                         self.colorticks = self.cbprop.colorticks
        559                         self.colorlabels = self.cbprop.colorlabels
    ~/plot_python/sample_unstructured/Plot_2D.py in __init__(self, arrays, colorticks, colorlabels, ranges, Nticks_list, max_find_method, tick_find_method, min_set, max_set, Ntick_set)
        826             elif max_find_method=='ceil':
        827                 tmp = maxval - np.array(ranges)
    --> 828                 plotmax = ranges[ np.where( tmp >= 0. )[0][-1]+1 ]
        829                 #print( 'ceil', plotmax )
    IndexError: index -1 is out of bounds for axis 0 with size 0
    1. Sorry about that. There have been several Plot_2D versions because it's continuously being developed. It looks like the default value of pretty_tick was set to True in your Plot_2D version. It happens when the values to be plotted are very low (<1e-6 or so) with the pretty_tick keyword turned on.

      You can try one of these:

      (1) Plot_2D( ds_fv['CO'][0,-1,:,:]*1e9 )   # convert CO mixing ratio: mol/mol --> ppbv

      (2) Plot_2D( ds_fv['CO'][0,-1,:,:], pretty_tick=False 

      Let me know if you have further questions/problems. 

      1. This resolved my issue, thank you!

  9. Anonymous


    I'm running into some issues with memory on Jupyter Hub when attempting to plot my MUSICAv0 outputs. I was successfully able to plot Ozone for June 2014, but when I use the same method to plot for July and August in that same Jupyter Notebook, I get both a File Save Error in JupyterHub and when I run – Plot_2D(ds_2['O3'][0,-1,:]*1e9, scrip_file=scrip_filename, pretty_tick=True, lon_range=[-130,-60], lat_range=[25,55], cmap=cm.hot_r, cmax=100, state=True, title='O3 Mixing Ratio over the US - July 2014', unit='ppbv', unit_offset=[0,-1]) – I get a memory error (see below). Thanks in advance for your help!

    MemoryError                               Traceback (most recent call last)
    <ipython-input-7-adb59e3aa7ba> in <module>
    ----> 1 Plot_2D(ds_2['O3'][0,-1,:]*1e9, scrip_file=scrip_filename, pretty_tick=True, lon_range=[-130,-60], lat_range=[25,55], cmap=cm.hot_r, cmax=100, state=True, title='O3 Mixing Ratio over the US - July 2014', unit='ppbv', unit_offset=[0,-1])
    ~/plot_python/f.e22.ne0CONUSne30x8_ne0CONUSne30x8_mt12.002_plot.v1/Plot_2D.py in __init__(self, var, lons, lats, lon_range, lat_range, scrip_file, ax, cmap, projection, grid_line, grid_line_lw, coast, country, state, resolution, feature_line_lw, feature_color, lonlat_line, lon_interval, lat_interval, font_family, label_size, colorbar, log_scale, log_scale_min, diff, orientation, shrink, pad, fraction, extend, colorticks, colorlabels, pretty_tick, nticks, cmax, cmin, title, title_size, title_bold, unit, unit_size, unit_bold, unit_italic, unit_offset, verbose)
        325                                                  np.array(self.lons_corners_add)), axis=0 )
        326             self.lats_corners = np.concatenate( (self.lats_corners, 
    --> 327                                                  np.array(self.lats_corners_add)), axis=0 )
        328             self.var = np.concatenate( (self.var, np.array(self.var_add)), axis=0 )
    <__array_function__ internals> in concatenate(*args, **kwargs)
    MemoryError: Unable to allocate array with shape (174369, 10, 1) and data type float64
    1. It looks like the required memory exceeded the available memory for plotting because the CONUS grid has a lot of grid cells and requires 10 lat/lon corners for each grid cell for plotting to keep original grid shape. I can think of two ways to solve this issue.
      (1) If you use Cheyenne/Casper on JupyterHub, specify the amount of memory as the maximum value (109 for Cheyenne and ~300? for Casper).
      (2) If you use other machines and you can't manage the memory, then I recommend you release the memory after you plot Ozone for June 2014. You may want to try "TMP = Plot_2D(...)" and "del TMP".

      Another way could be regridding the unstructured grid to the regular lat/lon structured grid for plotting. Please search "griddata" on this wiki page, it will show how to regrid unstructured grid to regular FV grid. In this way, you can save a lot of memory but lose the exact grid shape on the plot. But I believe generally it's fine for most of the applications to show scientific results.

      1. Thanks so much for your response Duseong!

        How would I go about specifying the amount of memory within JupyterHub? When I initially log onto JupyterHub it doesn't give me the option to launch a new job with adjusted settings. Is there a command I can use within JupyterHub to set the maximum limit?

        I appreciate your time.

        1. I am sorry that I wasn't aware of the new options on JupyterHub. It looks like you selected "login" node on the cluster selection page. Try to use "Cheyenne PBS batch" or "Casper PBS batch", then you will find this page below. 

          I recommend you use the values above I filled in. You need your project account. If you don't know your project account number, you can find it on https://sam.ucar.edu. Let me know if you have further questions!

          1. This worked for me really great. Thank you so much Duseong!