Basic Plotting#
BEFORE BEGINNING THIS EXERCISE - Check that your kernel (upper right corner, above) is NPL 2023b. This should be the default kernel, but if it is not, click on that button and select NPL 2023b.
This activity was developed primarily by Cecile Hannay and Jesse Nusbaumer.
For the atmospheric data, we will look at common variables in the atmospheric diagnostics. This notebook covers 3 basic plotting examples:
Exercise 1: Global lat/lon of surface temperature
Exercise 2: Zonal mean of short wave cloud forcing
Exercise 3: Temperature zonal mean with vertical levels
Some of the plotting in these examples are based on the AMWG Diagnostics Framework (ADF) and some are natively from the xarray functionality. xarray will be used for the data I/O, analysis, and some plotting, matplotlib and cartopy will aid in plotting, and numpy for calculations
import os
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cftime
import matplotlib as mpl
import matplotlib.path as mpath
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from matplotlib.gridspec import GridSpec
from matplotlib.lines import Line2D
from mpl_toolkits.axes_grid1 import make_axes_locatable
The first step is to grab an atmosphere (CAM) history file from your CESM model run
# Set your username here:
username = "PUT_USER_NAME_HERE"
# Here we point to the archive directory from your b1850.run_length simulation
monthly_output_path = f"/glade/derecho/scratch/{username}/archive/b1850.run_length/atm/hist"
# If you were unable to successfully run the b1850.run_length simulation, then feel free to use
# this provided simulation data instead:
#monthly_output_path = "/glade/campaign/cesm/tutorial/tutorial_2023_archive/b1850.run_length/atm/hist"
# Name of history file to plot
file_name = "b1850.run_length.cam.h0.0003-07.nc"
files = os.path.join(monthly_output_path, file_name)
files
'/glade/derecho/scratch/PUT_USER_NAME_HERE/archive/b1850.run_length/atm/hist/b1850.run_length.cam.h0.0003-07.nc'
ds = xr.open_dataset(files)
ds
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
File /glade/u/apps/opt/conda/envs/npl-2023b/lib/python3.10/site-packages/xarray/backends/file_manager.py:211, in CachingFileManager._acquire_with_cache_info(self, needs_lock)
210 try:
--> 211 file = self._cache[self._key]
212 except KeyError:
File /glade/u/apps/opt/conda/envs/npl-2023b/lib/python3.10/site-packages/xarray/backends/lru_cache.py:56, in LRUCache.__getitem__(self, key)
55 with self._lock:
---> 56 value = self._cache[key]
57 self._cache.move_to_end(key)
KeyError: [<class 'netCDF4._netCDF4.Dataset'>, ('/glade/derecho/scratch/PUT_USER_NAME_HERE/archive/b1850.run_length/atm/hist/b1850.run_length.cam.h0.0003-07.nc',), 'r', (('clobber', True), ('diskless', False), ('format', 'NETCDF4'), ('persist', False)), 'da467520-dd8b-4048-a691-3445d0f0f93e']
During handling of the above exception, another exception occurred:
FileNotFoundError Traceback (most recent call last)
Cell In[3], line 1
----> 1 ds = xr.open_dataset(files)
2 ds
File /glade/u/apps/opt/conda/envs/npl-2023b/lib/python3.10/site-packages/xarray/backends/api.py:566, in open_dataset(filename_or_obj, engine, chunks, cache, decode_cf, mask_and_scale, decode_times, decode_timedelta, use_cftime, concat_characters, decode_coords, drop_variables, inline_array, chunked_array_type, from_array_kwargs, backend_kwargs, **kwargs)
554 decoders = _resolve_decoders_kwargs(
555 decode_cf,
556 open_backend_dataset_parameters=backend.open_dataset_parameters,
(...)
562 decode_coords=decode_coords,
563 )
565 overwrite_encoded_chunks = kwargs.pop("overwrite_encoded_chunks", None)
--> 566 backend_ds = backend.open_dataset(
567 filename_or_obj,
568 drop_variables=drop_variables,
569 **decoders,
570 **kwargs,
571 )
572 ds = _dataset_from_backend_dataset(
573 backend_ds,
574 filename_or_obj,
(...)
584 **kwargs,
585 )
586 return ds
File /glade/u/apps/opt/conda/envs/npl-2023b/lib/python3.10/site-packages/xarray/backends/netCDF4_.py:590, in NetCDF4BackendEntrypoint.open_dataset(self, filename_or_obj, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta, group, mode, format, clobber, diskless, persist, lock, autoclose)
569 def open_dataset( # type: ignore[override] # allow LSP violation, not supporting **kwargs
570 self,
571 filename_or_obj: str | os.PathLike[Any] | BufferedIOBase | AbstractDataStore,
(...)
587 autoclose=False,
588 ) -> Dataset:
589 filename_or_obj = _normalize_path(filename_or_obj)
--> 590 store = NetCDF4DataStore.open(
591 filename_or_obj,
592 mode=mode,
593 format=format,
594 group=group,
595 clobber=clobber,
596 diskless=diskless,
597 persist=persist,
598 lock=lock,
599 autoclose=autoclose,
600 )
602 store_entrypoint = StoreBackendEntrypoint()
603 with close_on_error(store):
File /glade/u/apps/opt/conda/envs/npl-2023b/lib/python3.10/site-packages/xarray/backends/netCDF4_.py:391, in NetCDF4DataStore.open(cls, filename, mode, format, group, clobber, diskless, persist, lock, lock_maker, autoclose)
385 kwargs = dict(
386 clobber=clobber, diskless=diskless, persist=persist, format=format
387 )
388 manager = CachingFileManager(
389 netCDF4.Dataset, filename, mode=mode, kwargs=kwargs
390 )
--> 391 return cls(manager, group=group, mode=mode, lock=lock, autoclose=autoclose)
File /glade/u/apps/opt/conda/envs/npl-2023b/lib/python3.10/site-packages/xarray/backends/netCDF4_.py:338, in NetCDF4DataStore.__init__(self, manager, group, mode, lock, autoclose)
336 self._group = group
337 self._mode = mode
--> 338 self.format = self.ds.data_model
339 self._filename = self.ds.filepath()
340 self.is_remote = is_remote_uri(self._filename)
File /glade/u/apps/opt/conda/envs/npl-2023b/lib/python3.10/site-packages/xarray/backends/netCDF4_.py:400, in NetCDF4DataStore.ds(self)
398 @property
399 def ds(self):
--> 400 return self._acquire()
File /glade/u/apps/opt/conda/envs/npl-2023b/lib/python3.10/site-packages/xarray/backends/netCDF4_.py:394, in NetCDF4DataStore._acquire(self, needs_lock)
393 def _acquire(self, needs_lock=True):
--> 394 with self._manager.acquire_context(needs_lock) as root:
395 ds = _nc4_require_group(root, self._group, self._mode)
396 return ds
File /glade/u/apps/opt/conda/envs/npl-2023b/lib/python3.10/contextlib.py:135, in _GeneratorContextManager.__enter__(self)
133 del self.args, self.kwds, self.func
134 try:
--> 135 return next(self.gen)
136 except StopIteration:
137 raise RuntimeError("generator didn't yield") from None
File /glade/u/apps/opt/conda/envs/npl-2023b/lib/python3.10/site-packages/xarray/backends/file_manager.py:199, in CachingFileManager.acquire_context(self, needs_lock)
196 @contextlib.contextmanager
197 def acquire_context(self, needs_lock=True):
198 """Context manager for acquiring a file."""
--> 199 file, cached = self._acquire_with_cache_info(needs_lock)
200 try:
201 yield file
File /glade/u/apps/opt/conda/envs/npl-2023b/lib/python3.10/site-packages/xarray/backends/file_manager.py:217, in CachingFileManager._acquire_with_cache_info(self, needs_lock)
215 kwargs = kwargs.copy()
216 kwargs["mode"] = self._mode
--> 217 file = self._opener(*self._args, **kwargs)
218 if self._mode == "w":
219 # ensure file doesn't get overridden when opened again
220 self._mode = "a"
File src/netCDF4/_netCDF4.pyx:2464, in netCDF4._netCDF4.Dataset.__init__()
File src/netCDF4/_netCDF4.pyx:2027, in netCDF4._netCDF4._ensure_nc_success()
FileNotFoundError: [Errno 2] No such file or directory: '/glade/derecho/scratch/PUT_USER_NAME_HERE/archive/b1850.run_length/atm/hist/b1850.run_length.cam.h0.0003-07.nc'
Exercise 1: Make a lat-lon plot of TS#
To highlight plotting the variables from the CESM atmosphere (CAM) file, the first example will plot a simple global lat/lon plot of surface temperature TS
Grab data from first time stamp#
NOTE: This dataset has only one time
ts_0 = ds.TS.sel({"time": ds.TS.time.values[0]}).squeeze()
ts_0
The next step is to set up the map. Since we are plotting a global lat/lon, we will use the “Plate Carree” projection.
# define the colormap
cmap = mpl.colormaps["jet"]
# set up the figure with a Plate Carree projection
fig = plt.figure(figsize=(15, 10))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
# Plot the first timeslice of TS
img = ax.pcolormesh(ds.lon, ds.lat, ts_0, cmap=cmap, transform=ccrs.PlateCarree())
plt.title("Surface Temperature", fontsize=20)
# Set up colorbar
plt.colorbar(img, orientation="vertical", fraction=0.0242, pad=0.01)
plt.show()
Click here for the solution
Figure: Plotting solution.
Question:
The colorbar limits are set automatically by pcolormesh. How could you change the arguments for the pcolormesh function to set the plotting limits?
Click here for hints
Choose a maximum temperature of 300K and minimum temperature of 225K.
img = ax.pcolormesh(ds.lon, ds.lat, ts_0, vmax=300, vmin=225, cmap=cmap, transform=ccrs.PlateCarree())
Question:
How could we change the central longitude?
Click here for hints
The default central longitude is 0. Try setting it to 180. Then try other values from 0-360.ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree(central_longitude=180))
A second quick example is for xarray’s built-in plotting which uses the matplotlib and cartopy in the background. xarray makes creating a basic plot fairly simple.
# Xarray native plotting
# Set up figure and axis
fig, ax = plt.subplots(1, figsize=(20, 10))
# Plot the data straight from the xarray dataset
ts_0.plot.contourf(cmap="jet", levels=np.arange(220, 321, 5))
plt.show()
Click here for the solution
Figure: Plotting solution.
Exercise 2: Zonal plot of SWCF#
The second example will plot the short wave cloud forcing (SWCF) zonally.
Grab the variable data and mean over the the lon value
ds_swcf = ds.SWCF
# Get all the dataset dimensions
d = ds_swcf.dims
# Grab all dimensions to mean the lon values from
davgovr = [dim for dim in d if dim not in ("lev", "lat")]
# Make new dataset of zonal mean values
ds_swcf_zonal = ds_swcf.mean(dim=davgovr)
# print some values of new zonally-averaged SWCF variable
ds_swcf_zonal
# Create figure and axis
fig, ax = plt.subplots(1, figsize=(12, 7))
# Set Main title for subplots:
plt.title("Short Wave Cloud Forcing", fontsize=20)
# Plot value on y-axis and latitude on the x-axis
ax.plot(ds_swcf_zonal.lat, ds_swcf_zonal, c="goldenrod")
ax.set_xlim(
[max([ds_swcf_zonal.lat.min(), -90.0]), min([ds_swcf_zonal.lat.max(), 90.0])]
)
ax.set_xlabel("Latitude")
plt.show()
Click here for the solution
Figure: Plotting solution.
Question
What code could you add to set a legend for the plot line?
Click here for hints
enter the following lines after the ax.set_xlabel command and before the plt.show() command.
label = ds_swcf.time.values[0].strftime() # -> '0001-02-01 00:00:00'
line = Line2D([0], [0], label=label,
color="blue")
fig.legend(handles=[line],bbox_to_anchor=(-0.15, 0.15, 1.05, .102),loc="right",
borderaxespad=0.0,fontsize=16,frameon=False)
Question
What else could you label the line legend, if anything?
Click here for hints
Try setting the label value to something else like your simulation run name, the season or month, etc.
label = 'Season '
Question
How can you change the color on the plot and the legend?
Click here for hints
To change colors, try a different named color (e.g. “red” or “goldenrod”). You should make sure both the plot and the legend match or your plot won’t make sense.
Click here for a solution
Figure: Plotting solution.
Exercise 3: Plot of zonal T#
This example will plot the 3D zonal mean of temperature T with the pressure as the y-variable and latitude as the x-variable. We are showing you four different ways to make the same plot since each is valuable for plots with pressure on the y-axis and how to format those to get more information from the data.
# Remove all dimensions of length one
ds_t = ds.T.squeeze()
ds_t
Plot 1: Natively via xarray#
This plot uses the xarray native grids for the plot
y-axis is increasing with height
y-axis is not in log pressure
# Average over all dimensions except `lev` and `lat`.
d = ds_t.dims
davgovr = [dim for dim in d if dim not in ("lev", "lat")]
DS = ds_t.mean(dim=davgovr)
DS.transpose("lat", "lev", ...)
fig, ax = plt.subplots(1, figsize=(20, 10))
DS.plot.contourf(ax=ax, y="lev", cmap="Spectral_r")
plt.show()
Click here for a solution
Figure: Plotting solution.
Plot 2: Quick ADF style#
This plot uses the quick style of the AMWG Diagnostics Framework (ADF) packages.
y-axis is increasing with height
y-axis is not in log pressure
# Average over all dimensions except `lev` and `lat`.
d = ds_t.dims
davgovr = [dim for dim in d if dim not in ("lev", "lat")]
DS = ds_t.mean(dim=davgovr)
lev = DS["lev"]
lat = DS["lat"]
mlev, mlat = np.meshgrid(lev, lat)
# Generate zonal plot:
fig, ax = plt.subplots(1, figsize=(15, 7))
# Create zonal plot with vertical levels
img = ax.contourf(mlat, mlev, DS.transpose("lat", "lev"), cmap="Spectral_r")
# Format axis and ticks
ax.set_xlabel("Latitude")
fig.colorbar(img, ax=ax, location="right")
Click here for a solution
Figure: Plotting solution.
Questions
What differences do you see between the ADF quick plot and the xarray native plot?
Why might you want to use one or the other?
Where, vertically, do high pressures occur and what does this mean for where the ground would be located on this plot?
Click here for hints
Notice that the colorbar and axes labels are different between the two plots. Also note that the y-axis is going from lower to higher pressure values, but in the real atmosphere the highest pressures are near the surface (so the plot is inverted relative to the actual height of the layers above the Earth’s surface).
Plot 3: ADF style with reversed y-axis#
y-axis is decreasing with height
y-axis is not in log pressure
# Average over all dimensions except `lev` and `lat`.
d = ds_t.dims
davgovr = [dim for dim in d if dim not in ("lev", "lat")]
DS = ds_t.mean(dim=davgovr)
# print(DS.lev.min(),DS.lev.max())
lev = DS["lev"]
lat = DS["lat"]
mlev, mlat = np.meshgrid(lev, lat)
# Generate zonal plot:
fig, ax = plt.subplots(1, figsize=(15, 7))
# Create zonal plot with vertical levels
img = ax.contourf(mlat, mlev, DS.transpose("lat", "lev"), cmap="Spectral_r")
# Format axis and ticks
plt.gca().invert_yaxis()
ax.tick_params(which="minor", length=4, color="r")
# Set up colorbar
cbar_ax = fig.add_axes([0, 0, 0.1, 0.1])
posn = ax.get_position()
# Set position and size of colorbar position: [left, bottom, width, height]
cbar_ax.set_position([posn.x0 + posn.width + 0.005, posn.y0, 0.02, posn.height])
ax.set_xlabel("Latitude")
fig.colorbar(img, cax=cbar_ax, orientation="vertical")
Click here for a solution
Figure: Plotting solution.
Question:
What differences do you see on the y-axis between the result from Plot 2 and Plot 3?
Click here for hints
The plots are mirror images of each other because the y-axis has been flipped. Plot 2 has the highest pressures at the top of the y-axis while Plot 3 has the highest pressures at the bottom of the y-axis.
Question:
Where, vertically, do high pressures occur and what does this mean for where the ground would be located on these two figures ?
Click here for hints
The highest pressures in the atmosphere occur at the ground where the most atmospheric mass is pressing downward. This means that Plot 2 is oriented so the location where the ground is is the top of the plot, or up. In Plot 3 the orientation is so that the ground would be at the bottom of the plot, or down. Thus Plot 3 is often more intuitive to read because it’s oriented with the way we perceive the atmosphere.
Plot 4: Complex ADF style#
y-axis is decreasing with height
y-axis is in log pressure
# Average over all dimensions except `lev` and `lat`.
d = ds_t.dims
davgovr = [dim for dim in d if dim not in ("lev", "lat")]
DS = ds_t.mean(dim=davgovr)
# print(DS.lev.min(),DS.lev.max())
lev = DS["lev"]
lat = DS["lat"]
mlev, mlat = np.meshgrid(lev, lat)
# Generate zonal plot:
fig, ax = plt.subplots(1, figsize=(15, 7))
# Create zonal plot with vertical levels
img = ax.contourf(mlat, mlev, DS.transpose("lat", "lev"), cmap="Spectral_r")
# Format axis and ticks
plt.yscale("log")
plt.gca().invert_yaxis()
ax.tick_params(which="minor", length=4, color="r")
# Set up colorbar
cbar_ax = fig.add_axes([0, 0, 0.1, 0.1])
posn = ax.get_position()
# Set position and size of colorbar position: [left, bottom, width, height]
cbar_ax.set_position([posn.x0 + posn.width + 0.005, posn.y0, 0.02, posn.height])
ax.set_xlabel("Latitude")
fig.colorbar(img, cax=cbar_ax, orientation="vertical")
Click here for a solution
Figure: Plotting solution.
Question:
What differences do you see on the y-axis between the result from Plot 3 and Plot 4?
Click here for hints
The plots both have the y-axis oriented so that the highest pressures, or ground, is at the bottom of the plot. But the scale of the y-axis has changed so that in Plot 3 it is logarithmic.
Question:
What does that mean in terms of which part of the atmosphere are emphasized?
Click here for hints
With the logarithmic y-axis more of the upper atmosphere is shown on the figure. If you wanted to emphasize the boundary layer you might want to change the limits of the figure.