NCL_conOncon_2.py

This script illustrates the following concepts:

• Overlaying two sets of contours on a map

• Drawing the zero contour line thicker

• Changing the center longitude for a cylindrical equidistant projection

• Using a blue-white-red color map

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

• Original NCL script: https://www.ncl.ucar.edu/Applications/Scripts/conOncon_2.ncl

• Original NCL plot: https://www.ncl.ucar.edu/Applications/Images/conOncon_2_lg.png

Import packages:

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

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

Read in data:

# Open a netCDF data file using xarray default engine and load the data into xarrays
sst = xr.open_dataset(gdf.get("netcdf_files/sst8292a.nc"))
olr = xr.open_dataset(gdf.get("netcdf_files/olr7991a.nc"))

# Extract data for December 1982
sst = sst.isel(time=11, drop=True).SSTA
olr = olr.isel(time=47, drop=True).OLRA

# Fix the artifact of not-shown-data around 0 and 360-degree longitudes
sst = gv.xr_add_cyclic_longitudes(sst, 'lon')
olr = gv.xr_add_cyclic_longitudes(olr, 'lon')

Plot:

# Generate figure and axes
plt.figure(figsize=(8, 8))

# Set axes projection
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=-160))
ax.set_extent([100, 300, -60, 60], crs=ccrs.PlateCarree())

# Load in color map and specify contour levels
cmap = cmaps.BlWhRe
sst_levels = np.arange(-5.5, 6, 0.5)
# Draw SST contour
temp = sst.plot.contourf(ax=ax,
                         transform=ccrs.PlateCarree(),
                         cmap=cmap,
                         levels=sst_levels,
                         extend='neither',
                         add_colorbar=False,
                         add_labels=False,
                         zorder=0)
plt.colorbar(temp,
             ax=ax,
             orientation='vertical',
             ticks=np.arange(-5, 6, 1),
             drawedges=True,
             shrink=0.5,
             aspect=10)

# Draw map features on top of filled contour
ax.add_feature(cfeature.LAND, facecolor='lightgray', zorder=1)
ax.add_feature(cfeature.COASTLINE, edgecolor='gray', linewidth=0.5, zorder=1)

# Draw OLR contour
# Specify contour levels excluding 0
olr_levels = np.arange(-80, 0, 10)
olr_levels = np.append(olr_levels, np.arange(10, 50, 10))

rad = olr.plot.contour(ax=ax,
                       transform=ccrs.PlateCarree(),
                       levels=olr_levels,
                       colors='gray',
                       linewidths=0.5,
                       add_labels=False)
ax.clabel(rad, [-40, -20, 20], fmt='%d', inline=True, colors='black')

# Plot the zero contour with a black color
rad = olr.plot.contour(ax=ax,
                       transform=ccrs.PlateCarree(),
                       levels=[0],
                       colors='black',
                       linewidths=0.5,
                       add_labels=False)
ax.clabel(rad, [0], fmt='%d', inline=True, colors='black')

# Use geocat.viz.util convenience function to set axes tick values
gv.set_axes_limits_and_ticks(ax,
                             ylim=(-60, 60),
                             yticks=np.arange(-60, 90, 30),
                             xticks=np.arange(-80, 120, 30))

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

# Remove the degree symbol from tick labels
ax.yaxis.set_major_formatter(LatitudeFormatter(degree_symbol=''))
ax.xaxis.set_major_formatter(LongitudeFormatter(degree_symbol=''))

# Use geocat.viz.util convenience function to add minor and major tick lines
gv.add_major_minor_ticks(ax,
                         x_minor_per_major=3,
                         y_minor_per_major=3,
                         labelsize=10)

gv.set_titles_and_labels(ax,
                         maintitle=olr.long_name,
                         maintitlefontsize=14,
                         lefttitle='degC',
                         lefttitlefontsize=12,
                         righttitle='(W m s$^{-2}$)',
                         righttitlefontsize=12)
# Add center title
ax.text(0.35, 1.06, 'December 1982', fontsize=12, transform=ax.transAxes)

# Add lower text box
ax.text(1,
        -0.2,
        "CONTOUR FROM -80 TO 40 BY 10",
        horizontalalignment='right',
        transform=ax.transAxes,
        bbox=dict(boxstyle='square, pad=0.25',
                  facecolor='white',
                  edgecolor='black'))

plt.show()