NCL_sat_2.py

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 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")
/home/docs/checkouts/readthedocs.org/user_builds/geocat-examples/conda/latest/lib/python3.11/site-packages/xarray/coding/times.py:170: SerializationWarning: Ambiguous reference date string: 1-1-1 00:00:0.0. The first value is assumed to be the year hence will be padded with zeros to remove the ambiguity (the padded reference date string is: 0001-1-1 00:00:0.0). To remove this message, remove the ambiguity by padding your reference date strings with zeros.
  warnings.warn(warning_msg, SerializationWarning)

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 = 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()
mean Daily Sea Level Pressure, $\bf{SLP}$ $\bf{1963,}$ $\bf{January}$ $\bf{24th}$, hPa
/home/docs/checkouts/readthedocs.org/user_builds/geocat-examples/checkouts/latest/Gallery/MapProjections/NCL_sat_2.py:105: 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,

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

Gallery generated by Sphinx-Gallery