NCL_sat_2.py

This script illustrates the following concepts:

• Converting float data into short data

• Drawing filled contours over a satellite map

• Explicitly setting contour fill colors

• Finding local high pressure values

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

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

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

Import packages:

import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors
import matplotlib.ticker as mticker

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
ds = xr.open_dataset(gdf.get("netcdf_files/slp.1963.nc"), decode_times=False)

# Get data from the 21st timestep
pressure = ds.slp[21, :, :]

# Translate float values to short values
pressure = pressure.astype('float32')

# Convert Pa to hPa data
pressure = pressure * 0.01

# Fix the artifact of not-shown-data around 0 and 360-degree longitudes
wrap_pressure = gv.xr_add_cyclic_longitudes(pressure, "lon")

Create plot

# Set figure size
fig = plt.figure(figsize=(8, 8))

# Set global axes with a nearside perspective projection (equivalent to NCL's
# satellite projection)
proj = ccrs.NearsidePerspective(central_longitude=270.0,
                                central_latitude=45.0,
                                satellite_height=12742000)
ax = plt.axes(projection=proj)
ax.set_global()

# Add land, coastlines, and ocean features
ax.add_feature(cfeature.LAND, facecolor='lightgray', zorder=1)
ax.add_feature(cfeature.COASTLINE, linewidth=.3, zorder=2)
ax.add_feature(cfeature.OCEAN, facecolor='white')
ax.add_feature(cfeature.BORDERS, linewidth=.3)
ax.add_feature(cfeature.LAKES,
               facecolor='white',
               edgecolor='black',
               linewidth=.3)

# Create color map
colorvalues = [1020, 1036, 1500]
cmap = colors.ListedColormap(['None', 'lightgray', 'dimgrey'])
norm = colors.BoundaryNorm(colorvalues, 2)

# Plot contour data
p = wrap_pressure.plot.contourf(ax=ax,
                                zorder=2,
                                transform=ccrs.PlateCarree(),
                                levels=30,
                                cmap=cmap,
                                norm=norm,
                                add_labels=False,
                                add_colorbar=False)

p = wrap_pressure.plot.contour(ax=ax,
                               transform=ccrs.PlateCarree(),
                               linewidths=0.3,
                               levels=30,
                               cmap='black',
                               add_labels=False)

# low pressure contour levels- these will be plotted
# as a subscript to an 'L' symbol.
lowClevels = gv.find_local_extrema(pressure, lowVal=995, eType='Low')
highClevels = gv.find_local_extrema(pressure, highVal=1042, eType='High')

# Label regular contours with automatic matplotlib labeling
# Specify the levels to label every other contour level
ax.clabel(p,
          levels=np.arange(956, 1064, 8),
          inline=True,
          fontsize=12,
          colors='black',
          fmt="%.0f")

# Label low and high contours
gv.plot_extrema_labels(wrap_pressure,
                       ccrs.Geodetic(),
                       proj,
                       label_locations=lowClevels,
                       label='L')
gv.plot_extrema_labels(wrap_pressure,
                       ccrs.Geodetic(),
                       proj,
                       label_locations=highClevels,
                       label='H')

# Use gv function to set title and subtitles
gv.set_titles_and_labels(ax,
                         maintitle=r"$\bf{SLP}$" + " " + r"$\bf{1963,}$" + " " +
                         r"$\bf{January}$" + " " + r"$\bf{24th}$",
                         maintitlefontsize=20,
                         lefttitle="mean Daily Sea Level Pressure",
                         lefttitlefontsize=16,
                         righttitle="hPa",
                         righttitlefontsize=16)

# Set characteristics of text box
props = dict(facecolor='white', edgecolor='black', alpha=0.5)

# Place text box
ax.text(0.40,
        -0.1,
        'CONTOUR FROM 948 TO 1064 BY 4',
        transform=ax.transAxes,
        fontsize=16,
        bbox=props)

# Add gridlines to axis
gl = ax.gridlines(color='gray', linestyle='--')
gl.xlocator = mticker.FixedLocator(np.arange(-180, 180, 20))
gl.ylocator = mticker.FixedLocator(np.arange(-90, 90, 20))

# Make layout tight
plt.tight_layout()

plt.show()

/home/docs/checkouts/readthedocs.org/user_builds/geocat-examples/checkouts/latest/Gallery/MapProjections/NCL_sat_2.py:108: UserWarning: The following locations could not be translated into the desired projection: (100.0, -62.5). These locations will be dropped.
  gv.plot_extrema_labels(wrap_pressure,
/home/docs/checkouts/readthedocs.org/user_builds/geocat-examples/checkouts/latest/Gallery/MapProjections/NCL_sat_2.py:113: UserWarning: The following locations could not be translated into the desired projection: (100.0, 52.5). These locations will be dropped.
  gv.plot_extrema_labels(wrap_pressure,