# NCL_dev_2.py#

This script illustrates the following concepts:
• Calculating deviation from zonal mean

• Drawing zonal average plots

• Moving the contour informational label into the plot

• Changing the background color of the contour line labels

• Spanning part of a color map for contour fill

• Making the colorbar be vertical

• Paneling four subplots in a two by two grid using gridspec

• Changing the aspect ratio of a subplot

• Drawing color-filled contours over a cylindrical equidistant map

• Using a blue-white-red color map

See following URLs to see the reproduced NCL plot & script:

Import packages:

```import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LongitudeFormatter, LatitudeFormatter
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

import geocat.datafiles as gdf
import cmaps
import geocat.viz as gv
```

```# Open a netCDF data file using xarray default engine and load the data into xarrays
ds = xr.open_dataset(gdf.get("netcdf_files/b003_TS_200-299.nc"),
decode_times=False)

# Extract slice of data at first timestep
TS_0 = ds.TS.isel(time=0, drop=True)

# Calculate zonal mean
mean = TS_0.mean(dim='lon')

# Calculate deviation from time average by finding the temperatures averaged
# over all timesteps. Then that average is subtracted from the first timestep
time_avg = ds.TS.mean(dim='time')
time_dev = TS_0 - time_avg

# Fix the artifact of not-shown-data around 0 and 360-degree longitudes
```

Plot:

```# Specify projection for maps
proj = ccrs.PlateCarree()

# Generate figure (set its size (width, height) in inches)
fig = plt.figure(figsize=(8, 8))

# Create girdspec for layout, width_ratio is used to make the plots on the
# right narrower than the ones on the left
nrows=2,
width_ratios=[0.85, 0.15],
wspace=0.08)

# Create axis for plot with data from first timestep
ax1.coastlines(linewidths=0.25)

# Create axis for zonal mean temperature plot

# Create axis for deviation from time data plot
ax3.coastlines(linewidths=0.25)

# Create axis for colorbar

# Format ticks and ticklabels for the map axes
for ax in [ax1, ax3]:
# Use the geocat.viz function to set axes limits and ticks
gv.set_axes_limits_and_ticks(ax,
xlim=[-180, 180],
ylim=[-90, 90],
xticks=np.arange(-180, 181, 30),
yticks=np.arange(-90, 91, 30))

# Use the geocat.viz function to add minor ticks

# Use geocat.viz.util convenience function to make plots look like NCL
# plots by using latitude, longitude tick labels

# Removing degree symbol from tick labels to resemble NCL example
ax.yaxis.set_major_formatter(LatitudeFormatter(degree_symbol=''))
ax.xaxis.set_major_formatter(LongitudeFormatter(degree_symbol=''))

# Use the geocat.viz function to set axes limits and ticks for zonal average plot
gv.set_axes_limits_and_ticks(ax2,
xlim=[200, 310],
ylim=[-90, 90],
xticks=[200, 240, 280],
yticks=[])

# Use the geocat.viz function to add minor ticks to zonal average plot

# Plot contour lines for data at first timestep
contour = TS_0.plot.contour(ax=ax1,
transform=proj,
vmin=235,
vmax=305,
levels=np.arange(210, 311, 10),
colors='black',
linewidths=0.25,

# Label every other contour lines
ax1.clabel(contour, np.arange(220, 311, 20), fmt='%d', inline=True, fontsize=10)

# Set label backgrounds white
for txt in contour.labelTexts:

ax1.text(0.995,
0.03,
"CONTOUR FROM 210 TO 310 BY 10",
horizontalalignment='right',
transform=ax1.transAxes,
fontsize=8,
facecolor='white',
edgecolor='black'),
zorder=5)

# Add titles to top plot
size = 10
y = 1.05
ax1.set_title('Time(0)', fontsize=size, y=y)
ax1.set_title(TS_0.long_name, fontsize=size, loc='left', y=y)
ax1.set_title(TS_0.units, fontsize=size, loc='right', y=y)

# Plot zonal mean temperature
ax2.plot(mean.data, mean.lat, color='black', linewidth=0.5)

# Plot vertical reference line in zonal mean plot
ax2.axvline(273.15, color='black', linewidth=0.5)

# Import color map
cmap = cmaps.BlWhRe

# Truncate colormap to only use paler colors in the center of the colormap
cmap = gv.truncate_colormap(cmap, minval=0.22, maxval=0.74, n=14)

# Plot filled contour for deviation from time avg plot
deviations = time_dev.plot.contourf(ax=ax3,
transform=proj,
vmin=-14,
vmax=18,
levels=np.arange(-14, 20, 2),
cmap=cmap,

# Draw contour lines for deviation from time avg plot
time_dev.plot.contour(ax=ax3,
transform=proj,
vmin=-14,
vmax=18,
levels=np.arange(-14, 20, 2),
colors='black',
linewidths=0.25,
linestyles='solid',

# Add titles to bottom plot
ax3.set_title('Deviation from time ave', fontsize=size, y=y)
ax3.set_title(ds.TS.long_name, fontsize=size, loc='left', y=y)
ax3.set_title(ds.TS.units, fontsize=size, loc='right', y=y)

plt.colorbar(deviations,
cax=ax4,
ticks=np.linspace(-12, 16, 15),
drawedges=True)

plt.show()
```
```/home/docs/checkouts/readthedocs.org/user_builds/geocat-examples/conda/latest/lib/python3.11/site-packages/geocat/viz/util.py:730: UserWarning: Overwriting the cmap 'trunc(BlWhRe,0.22,0.74)' that was already in the registry.
mpl.colormaps.register(new_cmap, force=force)
```

Total running time of the script: (0 minutes 1.019 seconds)

Gallery generated by Sphinx-Gallery