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_vars(names=["time"])
ds = ds.isel(lev=0).drop_vars(names=["lev"])
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()