season plot with cartopy in python

1# Author: Daniel Rothenberg
2# Version: June 2, 2017
3
4import matplotlib.pyplot as plt
5plt.style.use(['seaborn-talk', 'seaborn-ticks'])
6
7import xbpch
8
9# First we read in a sample dataset containing GEOS-Chem output
10ds = xbpch.open_bpchdataset(
11    "/Users/daniel/workspace/bpch/test_data/ref_e2006_m2008.bpch",
12    diaginfo_file="/Users/daniel/Desktop/sample_nd49/diaginfo.dat",
13    tracerinfo_file="/Users/daniel/Desktop/sample_nd49/tracerinfo.dat",
14    dask=True, memmap=True
15)
16
17# Compute seasonal averages by doing splitting along the "seasons" corresponding
18# to each timestep, and taking the average over time in that group
19seasonal_o3 = (
20    ds['IJ_AVG_S_O3']
21    .isel(lev=0)  # select just surface values
22    .groupby('time.season').mean('time')
23)
24print(seasonal_o3)
25
26# Note that we now have a new dimension, "season", corresponding to the groups
27# we split the dataset by. We can use this dimension as a 'facet' to assemble
28# a collection of plots.
29# TODO: Cleanup axis proportions
30import cartopy.crs as ccrs
31g = seasonal_o3.plot.imshow(
32    x='lon', y='lat', # Use lat/lon for the x/y axis coordinates
33    vmin=0, vmax=60., cmap='gist_stern', # Colormap settings for all panels
34    col='season', col_wrap=2, # Facet over 'season', with 2 columns on the grid
35    transform=ccrs.PlateCarree(), # Geographic transform for coordinates
36    subplot_kws=dict(projection=ccrs.PlateCarree())
37        # Have each subpanel use this cartographic projection
38)
39for ax in g.axes.ravel():
40    ax.coastlines()
41plt.show()
42