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:
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import geocat.datafiles as gdf
import geocat.viz.util as gvutil
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import numpy as np

Import packages:

import xarray as xr
from matplotlib import colors

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 = gvutil.xr_add_cyclic_longitudes(pressure, "lon")

Create plot

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

# Set global axes with an orthographic projection
proj = ccrs.Orthographic(central_longitude=270, central_latitude=45)
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 = gvutil.findLocalExtrema(pressure, lowVal=995, eType='Low')
highClevels = gvutil.findLocalExtrema(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
gvutil.plotELabels(wrap_pressure,
                   ccrs.Geodetic(),
                   proj,
                   clabel_locations=lowClevels,
                   label='L')
gvutil.plotELabels(wrap_pressure,
                   ccrs.Geodetic(),
                   proj,
                   clabel_locations=highClevels,
                   label='H')

# Use gvutil function to set title and subtitles
gvutil.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()
mean Daily Sea Level Pressure, $\bf{SLP}$ $\bf{1963,}$ $\bf{January}$ $\bf{24th}$, hPa

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

Gallery generated by Sphinx-Gallery