NCL_lcmask_1.py

This script illustrates the following concepts:

• Drawing filled contours over a Lambert Conformal map

• Zooming in on a particular area on a Lambert Conformal map

• Creating a custom plot boundary

• Using a blue-white-red color map

• Setting contour levels using a min/max contour level and a spacing

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

• Original NCL script: http://ncl.ucar.edu/Applications/Scripts/lcmask_1.ncl

• Original NCL plot: http://ncl.ucar.edu/Applications/Images/lcmask_1_1_lg.png and http://ncl.ucar.edu/Applications/Images/lcmask_1_2_lg.png

Import packages:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
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 and disable time decoding due to missing necessary metadata
ds = xr.open_dataset(gdf.get("netcdf_files/atmos.nc"), decode_times=False)
# Extract a slice of the data
ds = ds.isel(time=0, drop=True)
ds = ds.isel(lev=0, drop=True)
V = ds.V
# Ensure longitudes range from 0 to 360 degrees
V = gv.xr_add_cyclic_longitudes(V, "lon")

Plot unmasked data:

# Generate figure and projection using Cartopy
plt.figure(figsize=(7, 10))
proj = ccrs.LambertConformal(central_longitude=0, standard_parallels=(45, 89))
# Set axis projection
ax = plt.axes(projection=proj, frameon=False)
# Set extent to include all longitudes and the northern hemisphere
ax.set_extent((0, 359, 0, 89), crs=ccrs.PlateCarree())
ax.coastlines(linewidth=0.5)

# Plot data and create colorbar
newcmp = cmaps.BlWhRe

wind = V.plot.contourf(
    ax=ax, cmap=newcmp, transform=ccrs.PlateCarree(), add_colorbar=False, levels=24
)
cbar = plt.colorbar(
    wind,
    ax=ax,
    orientation='horizontal',
    drawedges=True,
    ticks=np.arange(-48, 48, 8),
    pad=0.1,
    aspect=12,
)
cbar.ax.tick_params(length=0)  # remove tick marks but leave in labels

# Use geocat.viz.util convenience function to add left and right titles
gv.set_titles_and_labels(
    ax,
    lefttitle=V.long_name,
    lefttitlefontsize=16,
    righttitle=V.units,
    righttitlefontsize=16,
)

plt.show()

Mask data

masked = V.where(V.lat > 20)
masked = masked.where(masked.lat < 80)
masked = masked.where(masked.lon > 90)
masked = masked.where(masked.lon < 220)

# Rotate data to match NCL example
masked['lon'] = masked['lon'] + 180

Plot masked data

# Generate figure and projection using Cartopy
plt.figure(figsize=(10, 7))
proj = ccrs.LambertConformal(central_longitude=-22.5, standard_parallels=(45, 89))
# Set axis projection
ax = plt.axes(projection=proj)
ax.coastlines(linewidth=0.5)

# Make a custom boundary using convenience function
gv.set_map_boundary(ax, [-85, 40], [20, 80], south_pad=1)

# Plot data and create colorbar
wind = masked.plot.contourf(
    ax=ax, cmap=newcmp, transform=ccrs.PlateCarree(), add_colorbar=False, levels=24
)
cbar = plt.colorbar(
    wind,
    ax=ax,
    orientation='horizontal',
    drawedges=True,
    ticks=np.arange(-40, 44, 4),
    pad=0.1,
    aspect=18,
)
cbar.ax.tick_params(length=0)  # remove tick marks but leave in labels

# Use geocat.viz.util convenience function to add left and right titles
gv.set_titles_and_labels(
    ax,
    lefttitle=V.long_name,
    lefttitlefontsize=16,
    righttitle=V.units,
    righttitlefontsize=16,
)
plt.show()