metpy cross_section function gives nan in output

I'm trying to use metpy to take a cross-section of oceanographic data, because I'm not managing with the OceanSpy package and don't know if there are other packages more suitable for doing this. Advice is welcome :-)

I downloaded some sample data from CMEMS manually selecting a small map area and ticking the variables uo, vo and zos from https://data.marine.copernicus.eu/product/GLOBAL_MULTIYEAR_PHY_001_030/download?dataset=cmems_mod_glo_phy_my_0.083_P1M-m

I read the data using xarray and parse through metpy, focussing on a single timestep for simplification:

ds = xr.open_mfdataset('/Users/0448257/Data/cmems_mod_glo_phy-cur_anfc_0.083deg_P1M.nc')
data = ds.isel(time=30).metpy.parse_cf().squeeze()
print(data)

<xarray.Dataset>
Dimensions:    (depth: 35, latitude: 174, longitude: 175)
Coordinates:
  * depth      (depth) float32 0.494 1.541 2.646 3.819 ... 643.6 763.3 902.3
  * latitude   (latitude) float32 -67.33 -67.25 -67.17 ... -53.08 -53.0 -52.92
    time       datetime64[ns] 2023-05-16T12:00:00
  * longitude  (longitude) float32 -72.0 -71.92 -71.83 ... -57.67 -57.58 -57.5
    metpy_crs  object Projection: latitude_longitude
Data variables:
    vo         (depth, latitude, longitude) float32 dask.array<chunksize=(35, 174, 175), meta=np.ndarray>
    uo         (depth, latitude, longitude) float32 dask.array<chunksize=(35, 174, 175), meta=np.ndarray>

Here's a plot of where I'm trying to do the cross-section:

start = (-68, -56) #(-66, -60) 
end = (-64, -64) #(-63, -62) 
fig, ax = plt.subplots()
ax.contourf(data['longitude'], data['latitude'], data['uo'][0,:,:])
ax.plot((start[0],end[0]), (start[1], end[1]), marker='o', color='black')

Plot of u data with cross-section

I calculate the cross-section without errors, but when I print or plot the data there are nan in it...

It makes sense to get some nans deeper down along the coasts as the ocean narrows there, but not in the middle of the cross-section. I've also tried taking a smaller cross-section to make sure I never have NA in the data, but that didn't help. Setting the interpolation type to nearest also didn't help.

Any ideas?

Some sample output:

cross = cross_section(data, start, end).set_coords(('latitude', 'longitude'))
print(cross)

<xarray.Dataset>
Dimensions:    (depth: 35, index: 100)
Coordinates:
  * depth      (depth) float32 0.494 1.541 2.646 3.819 ... 643.6 763.3 902.3
    time       datetime64[ns] 2023-05-16T12:00:00
    metpy_crs  object Projection: latitude_longitude
    longitude  (index) float64 -56.0 -56.09 -56.19 ... -63.86 -63.93 -64.0
    latitude   (index) float64 -68.0 -67.96 -67.92 ... -64.08 -64.04 -64.0
  * index      (index) int64 0 1 2 3 4 5 6 7 8 9 ... 91 92 93 94 95 96 97 98 99
Data variables:
    vo         (depth, index) float32 dask.array<chunksize=(35, 100), meta=np.ndarray>
    uo         (depth, index) float32 dask.array<chunksize=(35, 100), meta=np.ndarray>

print(cross.uo.values)

[[        nan         nan         nan ... -0.00716022  0.01993629
   0.04136517]
 [        nan         nan         nan ... -0.00686111  0.02008043
   0.04135653]
 [        nan         nan         nan ... -0.00662968  0.02014319
   0.04128024]
 ...

This is a case of the horrible cross between normal points being (x, y), but geographically using (latitude, longitude). In MetPy for a lot of cases, we try to adhere (at least for individual function arguments) to longitude, then latitude. HOWEVER, we stuck with (latitude, longitude) points for cross_section. I'm not sure how best to deal with this problem in a way that avoids user problems.

At any rate, if I change your code to flip the order of the entries in start and end:

cross = cross_section(data, start[::-1], end[::-1])

I only get nan at the bottom

array([[0.286874  , 0.29746435, 0.27922769, ..., 0.06116014, 0.03308784,
        0.0225837 ],
       [0.28138065, 0.29160811, 0.27359135, ..., 0.05810829, 0.03002688,
        0.01892148],
       [0.27649769, 0.28672514, 0.26870839, ..., 0.05510405, 0.0275854 ,
        0.01648   ],
       ...,
       [       nan,        nan,        nan, ...,        nan,        nan,
               nan],
       [       nan,        nan,        nan, ...,        nan,        nan,
               nan],
       [       nan,        nan,        nan, ...,        nan,        nan,
               nan]])

and the plot looks reasonable (plt.imshow(cross.uo.values)):