Converting lat/lon coordinates into meters/kilometers for use in a metpy HRRR cross section

I am trying to create custom cross sections of archived HRRR Grib2 output data. I had been following the cross section example provided here and followed up on all issues I had with the file format itself also on the unidata site here. I have produced some cross section plots such as the following where my x-axis utilizes latitude and y-axis utilizes isobaric pressure as seen in the plot here:

My goal is to output my plots with the x-axis showing distance from the start of my transect line to the end of the transect line. This would help me determine the horizontal scale of certain near-surface meteorological features including lakebreeze, outflow, etc. An example of what I am hoping to do is in the following photo, where the x axis indicates distance along the transect line in meters or km instead of gps coordinates:

How can I convert the coordinates to distances for my plot?

My Code:

#input points for the transect line via longitude/latitude coordinates
startpoint = (42.857, -85.381)
endpoint = (43.907, -83.910)

# Import
grib = pygrib.open('file.grib2')

#use grib message to apply lat/long info to data set
msg = grib.message(1) 

#Convert grib file into xarray and assign x and y coordinate values
#(HRRR utilizes lambert_conformal_conic by default remember this)
ds = xr.open_dataset(file, engine="cfgrib",
                     backend_kwargs={'filter_by_keys':{'typeOfLevel': 'isobaricInhPa'}})

ds = ds.metpy.assign_crs(CRS(msg.projparams).to_cf()).metpy.assign_y_x()

#metpy cross section function to create a transect line and create cross section.
cross = cross_section(ds, startpoint, endpoint).set_coords(('latitude', 'longitude'))
 
#create variables
temperature = cross['t']
pressure = cross['isobaricInhPa']
cross['Potential_temperature'] = mpcalc.potential_temperature(cross['isobaricInhPa'],cross['t'])
cross['u_wind'] = cross['u'].metpy.convert_units('knots')
cross['v_wind'] = cross['v'].metpy.convert_units('knots')
cross['t_wind'], cross['n_wind'] = mpcalc.cross_section_components(cross['u_wind'],cross['v_wind'])
cross['qv'] = cross['q'] *1000* units('g/kg')

#HRRR plot test
fig = plt.figure(1, figsize=(20,9))
ax = plt.axes()
#levels = np.linspace(325,365,50)
temp = ax.contourf(cross['latitude'], cross['isobaricInhPa'], cross['qv'], 100, cmap='rainbow')
clb = fig.colorbar(temp)
clb.set_label('g $\mathregular{kg^{-1}}$')
theta_contour = ax.contour(cross['latitude'], cross['isobaricInhPa'], cross['Potential_temperature'],
                           400, colors='k', linewidths=2)
theta_contour.clabel(theta_contour.levels[1::2], fontsize=8, colors='k', inline=1,
                     inline_spacing=8, fmt='%i', rightside_up=True, use_clabeltext=True)
ax.set_ylim(775,1000)
ax.invert_yaxis()
plt.title('HRRR contour fill of Mixing ratio(g/kg), contour of Potential Temperature (K),\n Tangential/Normal Winds (knots)')
plt.title('Run: '+date+'', loc='left', fontsize='small')
plt.title('Valid: '+date+' '+f_hr, loc='right', fontsize='small')
plt.xlabel('Latitude')
plt.ylabel('Pressure (hPa)')
wind_slc_vert = list(range(0, 19, 2)) + list(range(19, 29))
wind_slc_horz = slice(5, 100, 5)
ax.barbs(cross['latitude'][wind_slc_horz], cross['isobaricInhPa'][wind_slc_vert],
         cross['t_wind'][wind_slc_vert, wind_slc_horz],
         cross['n_wind'][wind_slc_vert, wind_slc_horz], color='k')
 
# Adjust y-axis to log scale
ax.set_yscale('symlog')
ax.set_yticklabels(np.arange(1000, 775,-100))
#ax.set_ylim(cross['isobaricInhPa'].max(), cross['isobaricInhPa'].min())
ax.set_yticks(np.arange(1000, 775,-100))

plt.show()

You should be able to do this using PyPROJ's Geod class, which is what MetPy uses under the hood to calculate e.g. lat_lon_grid_deltas. Something like this I think will work:

import pyproj
geod = pyproj.Geod(ellps='sphere')
_, _, dist = geod.inv(cross['longitude'][0], cross['latitude'][0],
                      cross['longitude'], cross['latitude'])

I haven't actually tried that code, so you may need to convert the xarray DataArrays into numpy plain arrays.