How do I overlay a line (e.g. Shapely LineString) onto a plot?

[ocoudert]

Hello,

I am trying to overlay a line onto a plot, for example a segment of the trajectory of this year's solar eclipse on contour plots.

The first attempt was with PlotGeometry, but I don't understand how the points of the line get geocoded, or how the CRS is specified.

gLine = geodL.inv_intermediate(start[1], start[0], end[1], end[0], step, return_back_azimuth=True)
lineShape = LineString(list(zip(gLine.lons, gLine.lats)))
linePlt = metpy.plots.PlotGeometry()
linePlt.geometry = lineShape

This fails with traitlets.traitlets.TraitError: The 'geometry' trait of a PlotGeometry instance expected an Iterable, not the LineString <LINESTRING (...

The second attempt was using PlotVector, but I couldn't get this to work:

lineValue = [0 for x in range(len(gLine.lons))] # filler data, for debug
lineStationIDs = [x for x in range(len(gLine.lons))] # filler data, for debug
lineDF = pd.DataFrame({ 'stid' : lineStationIDs, 'data' : lineValue, 'longitude' : gLine.lons, 'latitude' : gLine.lats })
lineDA = lineDF.to_xarray()
lineDA = lineDA.assign_attrs({ 'time' : dsTime })
lineDA = lineDA.assign_attrs({ 'name' : 'Houston to Cleveland' })
lineDA = lineDA.set_coords(['longitude','latitude'])
lineDA['longitude'] = units.degrees_east * lineDA['longitude']
lineDA['latitude'] = units.degrees_north * lineDA['latitude']
lineDA = lineDA.metpy.assign_crs(
    grid_mapping_name = proj
)
linePlt = metpy.plots.PlotVector() 
linePlt.data = lineDA
linePlt.color = 'y'
linePlt.fill = 'y'
linePlt.marker = '.-'

This results with AttributeError: 'PlotVector' object has no attribute '_build'

Do I need to merge all of the desired data into a single dataset? If so, how?

Setup:

• Python 3.11.4
• MetPy 1.6.1
• ecCodes 2.27.0
  These are set by company IT, I cannot upgrade.

[ocoudert]

Edited to include the line
linePlt.geometry = lineShape
in the first attempt. This is the line causing the error.

[ocoudert]

Here is the full code for context:

import xarray as xr
import pandas as pd
import metpy
from metpy.units import units
import matplotlib.pyplot as plt 
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from metpy.plots.declarative import (BarbPlot, ContourPlot,
                                     FilledContourPlot, MapPanel,
                                     PanelContainer, PlotObs)
from metpy.interpolate import cross_section
from shapely.geometry import LineString

# Create the dataset
filePath = 'C:\\temp\\gfs.t00z.pgrb2.0p25.f004'
ds = xr.open_dataset(filePath, engine='cfgrib',filter_by_keys={'typeOfLevel': 'isobaricInhPa'})
ds = ds.metpy.parse_cf().squeeze()
ds = ds.metpy.sel(latitude=slice(70, 10), longitude=slice(360 - 150, 360 - 55)) # focus on CONUS only
dtime=ds['time'].values
dsTime = dtime.astype('datetime64[us]').tolist() # convert to Python datetime

# Set units
ds['t'] = units('kelvin')*ds['t']
ds['r'] = units('percent')*ds['r']
ds['gh'] = units('m^2/s^2')*ds['gh']

# Select param to plot
param = 't'
paramRange = (200,330,20)
paramUnits = str(ds[param].metpy.units)
levelhPa = 900
colorMap = 'BuPu'

# Create a line (non-CF data)
start = (29.762778, -95.383056) # Houston
end = (41.482222, -81.669722) # Cleveland
# get projection and get points along line between start and end
if ds[param].metpy.cartopy_crs.name == 'unknown':
    proj = ccrs.PlateCarree()
else:
    proj = ds[param].metpy.cartopy_crs
geodL = proj.get_geod()
step = 20
gLine = geodL.inv_intermediate(start[1], start[0], end[1], end[0], step, return_back_azimuth=True)
# Verify path with Matplotlib
axL = plt.axes(projection=proj)
axL.set_extent([-125, -74,20, 55], proj)
axL.coastlines()
axL.plot(gLine.lons, gLine.lats, '.-', color='yellow', linewidth=1.5, markersize=3)
plt.show() # it works!
plt.close() 

# Attempt 1
# using Shapely LineString
lineShape = LineString(list(zip(gLine.lons, gLine.lats)))
linePlt = metpy.plots.PlotGeometry()
# How is the data geocoded?  How is CRS specified?
# Following line fails with traitlets.traitlets.TraitError: The 'geometry' trait of a PlotGeometry instance expected an Iterable, not the LineString <LINESTRING (...
##linePlt.geometry = lineShape

# Attempt 2
# Create line dataset - easiest way is to create a DataFrame and convert it to a DataArray
lineValue = [0 for x in range(len(gLine.lons))] # filler data, for debug
lineStationIDs = [x for x in range(len(gLine.lons))] # filler data, for debug
lineDF = pd.DataFrame({ 'stid' : lineStationIDs, 'data' : lineValue, 'longitude' : gLine.lons, 'latitude' : gLine.lats })
lineDA = lineDF.to_xarray()
lineDA = lineDA.assign_attrs({ 'time' : dsTime })
lineDA = lineDA.assign_attrs({ 'name' : 'Houston to Cleveland' })
lineDA = lineDA.set_coords(['longitude','latitude'])
lineDA['longitude'] = units.degrees_east * lineDA['longitude']
lineDA['latitude'] = units.degrees_north * lineDA['latitude']
lineDA = lineDA.metpy.assign_crs(
    grid_mapping_name = proj
)
linePlt = metpy.plots.PlotVector() # end result: AttributeError: 'PlotVector' object has no attribute '_build'
linePlt.data = lineDA
linePlt.color = 'y'
linePlt.fill = 'y'
linePlt.marker = '.-'

# Define contours
cntr2 = ContourPlot()
cntr2.data = ds
cntr2.field = param
cntr2.level = levelhPa * units.hPa
cntr2.contours = list(range(paramRange[0],paramRange[1],paramRange[2]))
cntr2.linecolor = 'black'
cntr2.linestyle = 'solid'
cntr2.clabels = True

# Fill contours
cfill = FilledContourPlot()
cfill.data = ds
cfill.field = param
cfill.level = levelhPa * units.hPa
cfill.contours = list(range(paramRange[0],paramRange[1],paramRange[2]))
cfill.colormap = colorMap
cfill.colorbar = 'horizontal'
cfill.plot_units = paramUnits

# Create panel for CONUS
panel = MapPanel()
panel.area = [-125, -74, 20, 55] # CONUS
panel.projection = 'lcc'
panel.layers = ['states', 'coastline', 'borders']
panel.title = f'{cfill.level.m} hPa {param} ({paramUnits}) at {dsTime}'
print('DEBUG: got this far...')
panel.plots = [cfill, cntr2, linePlt] 
pc = PanelContainer()
pc.size = (15,15)
pc.panels = [panel]

pc.show()

[ocoudert]

The error in the code above is that geometry should be specified as a list. The geometry attribute is expecting a collection of shapes.
The fix is simple:

linePlt.geometry = [lineShape]