NCL_conwomap_5.py

This script illustrates the following concepts:

• Drawing a simple contour plot

• Making an axis logarithmic in a contour plot

• Changing the labels and tickmarks on a contour plot

• Creating a main title

• Using the geocat-comp method interp_hybrid_to_pressure

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

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

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

Import packages:

import matplotlib.pyplot as plt
from matplotlib.ticker import ScalarFormatter
import numpy as np
import xarray as xr
import cmaps

from geocat.comp import interp_hybrid_to_pressure
import geocat.viz as gv
import geocat.datafiles as gdf

Read in data:

# Open a netCDF data file using xarray default engine and load the data into xarrays

ds = xr.open_dataset(gdf.get("netcdf_files/atmos.nc"), decode_times=False)

# Extract the data needed
u = ds.U[0, :, :, :]  # U component of wind
hyam = ds.hyam  # hybrid A coefficient
hybm = ds.hybm  # hybrid B coefficient
ps = ds.PS  # surface pressures in Pascals
p0 = 100000  # surface reference pressure in Pascals

# Specify output pressure levels
new_levels = np.array([1000, 950, 800, 700, 600, 500, 400, 300,
                       200])  # in millibars
new_levels = new_levels * 100  # convert to Pascals

# Interpolate pressure coordinates form hybrid sigma coord
u_int = interp_hybrid_to_pressure(u,
                                  ps[0, :, :],
                                  hyam,
                                  hybm,
                                  p0=p0,
                                  new_levels=new_levels,
                                  method='log')

# Calculate zonal mean of u component of wind
uzon = u_int.mean(dim='lon')

Plot:

# Generate figure (set its size (width, height) in inches)
plt.figure(figsize=(7, 10))
ax = plt.axes()

# Format log axis
plt.yscale('log')
ax.yaxis.set_major_formatter(ScalarFormatter())

# Use geocat.viz.util convenience function to set axes parameters
gv.set_axes_limits_and_ticks(ax,
                             ylim=(20000, 100000),
                             yticks=[100000, 70000, 50000, 30000],
                             yticklabels=['1000', '700', '500', '300'],
                             xticks=np.arange(-60, 90, 30),
                             xticklabels=['60S', '30S', '0', '30N', '60N'])

# Us geocat.viz.util convenience function to add minor and major ticks
gv.add_major_minor_ticks(ax,
                         x_minor_per_major=3,
                         y_minor_per_major=0,
                         labelsize=16)

# Specify colormap
newcmap = cmaps.ncl_default

# Plot filed contours
p = uzon.plot.contourf(ax=ax,
                       levels=13,
                       vmin=-8,
                       vmax=40,
                       cmap=newcmap,
                       add_colorbar=False,
                       add_labels=False)

# Plot contour lines
uzon.plot.contour(ax=ax,
                  levels=13,
                  vmin=-8,
                  vmax=40,
                  colors='black',
                  linewidths=0.5,
                  linestyles='solid',
                  add_labels=False)

# Create colorbar
cbar = plt.colorbar(p,
                    ax=ax,
                    drawedges=True,
                    extendrect=True,
                    extendfrac='auto',
                    ticks=np.arange(-8, 44, 4),
                    orientation='horizontal',
                    pad=0.075,
                    aspect=11)

# Set colorbar tick label size
cbar.ax.tick_params(labelsize=14)

# Use geocat.vix convenience function to set titles and labels
gv.set_titles_and_labels(ax,
                         maintitle="Logarithmic axis",
                         maintitlefontsize=18,
                         lefttitle="Zonal Wind",
                         lefttitlefontsize=16)

# Show plot
plt.show()

/home/docs/checkouts/readthedocs.org/user_builds/geocat-examples/conda/latest/lib/python3.11/site-packages/geocat/comp/interpolation.py:133: UserWarning: Interpolation point out of data bounds encountered
  return func_interpolate(new_levels, xcoords, data, axis=interp_axis)