Interpolation breaks (upgrade from 0.17.0 to 0.19.0)

I am using xarray to sample a 3D field along 2D Lagrangian trajectories (output from Parcels). The trajectory data consists of several spatial variables (longitude, latitude, depth), stored along the dimensions obs (essentially a time index) and traj (trajectory index). My aim is to interpolate the 3D field along 2D Lagrangian trajectories and obtain along-trajectory depth profiles of the data.

The code I used:

# Obtain profiles along all trajectories at the 0th observation
# `EulerianField` has dimensions ('Z', 'YC', 'XC')
# `LagrangianTrajectories.lon.isel(obs=0)` has dimension (`traj`); same for `lat`
EulerianField.interp(XC=LagrangianTrajectories.isel(obs=0).lon, 
                     YC=LagrangianTrajectories.isel(obs=0).lat, 
                     method='linear')

What happens in xarray 0.17.0:
Everything runs fine and quickly. A dataArray of the field is returned with dimensions (Z, traj); a depth profile for each trajectory (at obs=0).

What happens in xarray 0.19.0:
I get the following error:

numpy.core._exceptions.MemoryError: Unable to allocate 179. GiB for an array with shape (49, 62590, 62590) and data type bool

Note that the length of Z is 49, while that of traj is 62590. I suspect that xarray somehow tries to return an array with dimensions (Z, YC, XC), with YC and XC now having length 62590.

Minimal Complete Verifiable Example:
I haven't been able to reduce the issue to an MCVE. When I subset the trajectory dataset using .isel(obs=0, traj=slice(0,10), the behavior is as expected again (returning a dataArray with dimensions (Z, traj). I can send the files I'm working with using WeTransfer for testing, but they're several GBs. I hope that the issue contains enough hints for determining the cause.

Anything else we need to know?:
The trajectory data contains some NaNs, but when there are NaNs in the subsetted (traj=slice(0,10)), everything still works fine, so I don't think that's the cause.

Environment:

Start at the interpolation call, comparing behavior between xarray 0.17.0 and 0.19.0 with the reported dimensions: Z=49 and traj=62590. Reduce the trajectory input or obtain the referenced files to reproduce the memory allocation, then verify that interpolation returns dimensions (Z, traj) without attempting the 179 GiB array.