Plotting lat and lon on satellite image using MetPy's Declarative syntax

Question

I am attempting to find a way to visualize the separate regions/phases of the MJO. I believe one way to do so would be by plotting the longitude lines that separate each phase region (at roughly 60E, 80E, 100E, 120E, 140E, 160E, 180), but I am unsure if it is possible to add to my existing plots.

I am using GRID-Sat B1 data from NCEI. Here is what my current code looks like:

import matplotlib.pyplot as plt
from metpy.plots import declarative, colortables
import cartopy.crs as ccrs
import xarray as xr

file = "GRIDSAT-B1.2003.11.23.00.v02r01.nc"
dataset = xr.open_dataset(file)
vtime = dataset.time.values.astype('datetime64[s]').astype('O')
date_long = vtime[0]
date = date_long.strftime("%d-%b-%Y-%HZ")

# Create water vapor image
img = declarative.ImagePlot()
img.data = dataset
img.field = 'irwvp'
img.colormap = 'WVCIMSS_r'
img.image_range = (180, 280)
panel = declarative.MapPanel()
panel.layers = ['coastline', 'borders']
panel.title = f'GridSat-B1 (Water Vapor Imagery): {date}'
panel.projection = (ccrs.Mollweide(central_longitude=-240))
panel.area = ([-370, -140, -30, 30])
panel.layout = (2, 1, 2)
panel.plots = [img]

# Create the IR image
img2 = declarative.ImagePlot()
img2.data = dataset
img2.field = 'irwin_cdr'
img2.colormap = 'turbo_r' #maybe use cubehelix instead?
img2.image_range = (180, 300)
panel2 = declarative.MapPanel()
panel2.layers = ['coastline', 'borders']
panel2.title = f'GridSat-B1 (Infrared Imagery): {date}'
panel2.projection = (ccrs.Mollweide(central_longitude=-240))
panel2.area = ([-370, -140, -30, 30])
panel2.layout = (2, 1, 1)
panel2.plots = [img2]

# Plot both panels in one figure
pc = declarative.PanelContainer()
pc.size = (20, 14)
pc.panels = [panel, panel2]
pc.show()

Here is the current output that is created when I run the script: Nov03.png

Any help/suggestions are appreciated - thanks in advance!

Answer 1

There's nothing built into MetPy's declarative interface, but fortunately the MapPanel objects expose a .ax attribute that gives you a Matplotlib Axes object and all its plotting methods:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import metpy.plots as mpplots
import numpy as np
import xarray as xr

file = "/Users/rmay/Downloads/GRIDSAT-B1.2003.11.23.00.v02r01.nc"
dataset = xr.open_dataset(file)
vtime = dataset.time.values.astype('datetime64[s]').astype('O')
date_long = vtime[0]
date = date_long.strftime("%d-%b-%Y-%HZ")

# Create water vapor image
img = mpplots.ImagePlot()
img.data = dataset
img.field = 'irwvp'
img.colormap = 'WVCIMSS_r'
img.image_range = (180, 280)
panel = mpplots.MapPanel()
panel.layers = ['coastline', 'borders']
panel.title = f'GridSat-B1 (Water Vapor Imagery): {date}'
panel.projection = ccrs.Mollweide(central_longitude=-240)
panel.area = (-370, -140, -30, 30)
panel.layout = (2, 1, 2)
panel.plots = [img]

# Create the IR image
img2 = mpplots.ImagePlot()
img2.data = dataset
img2.field = 'irwin_cdr'
img2.colormap = 'turbo_r' #maybe use cubehelix instead?
img2.image_range = (180, 300)
panel2 = mpplots.MapPanel()
panel2.layers = ['coastline', 'borders']
panel2.title = f'GridSat-B1 (Infrared Imagery): {date}'
panel2.projection = ccrs.Mollweide(central_longitude=-240)
panel2.area = (-370, -140, -30, 30)
panel2.layout = (2, 1, 1)
panel2.plots = [img2]

# Plot both panels in one figure
pc = mpplots.PanelContainer()
pc.size = (20, 14)
pc.panels = [panel, panel2]

lons = np.array([60, 80, 100, 120, 140, 160, 180]).reshape(1, -1)
lats = np.linspace(-90, 90).reshape(-1, 1)

# Match up the arrays into 2xN arrays fit to plot in call
lons, lats = np.broadcast_arrays(lons, lats)

# Needs to be *after* the panels are assigned to a PanelContainer
# Using Geodetic gives lines interpolated on the curved globe
panel.ax.plot(lons, lats, transform=ccrs.Geodetic(), color='black', linewidth=3)
panel2.ax.plot(lons, lats, transform=ccrs.Geodetic(), color='black', linewidth=3)

pc.show()

(Note: it's not recommended to import from metpy's declarative module directly since that's an implementation detail subject to change--just get things from metpy.plots). So this is using Matplotlib's standard call to plot to draw the lines. Another option would be to use CartoPy's Gridliner.