Trouble using advanced indexing through space and time via opendap?

[najascutellatus]

Hello, I'm trying to calculate metrics at each point along a glider track and compare to a few different ocean models at the nearest grid location. I'm attempting to use advanced indexing to find the nearest neighbor to each surfacing, even though its such a small dataset (1100 that I am trying to work with in the full track), it's painfully slow when using OpenDAP and will usually not complete due to timeouts. Am I doing something wrong or is the server slow? Is there a more efficient way of accessing these specific points? Would looping through the points and using .sel be a better option for this specific case?

For brevity, I downsampled the glider to just 10 surfacing. It's still very slow.

import xarray as xr
import datetime as dt

gtime = [dt.datetime(2021,6,13,6,0,0),
         dt.datetime(2021,6,23,6,0,0),
         dt.datetime(2021,7,3,6,0,0),
         dt.datetime(2021,7,13,6,0,0),
         dt.datetime(2021,7,23,6,0,0),
         dt.datetime(2021,8,2,6,0,0),
         dt.datetime(2021,8,12,6,0,0),
         dt.datetime(2021,8,22,6,0,0),
         dt.datetime(2021,9,1,6,0,0),
         dt.datetime(2021,9,11,6,0,0),
         dt.datetime(2021,9,21,6,0,0)
        ]
glon = [
    265.41368532029145,
    266.2228786882868, 
    267.2704577270577, 
    268.85892911601684, 
    269.18903847079787, 
    269.3574762286588, 
    270.30725702433676, 
    271.4906198819087, 
    272.21700585406757, 
    272.0912898757826, 
    271.76470859324746
]
glat = [
    27.638446297249928, 
    27.056615139362126, 
    27.15419373137167,
    27.40356144679907,
    27.43763012714433,
    26.666685239752848,
    27.774503201623297,
    28.2032588731131,
    28.11229420537534,
    28.651781244954552, 
    29.194729088837246
]

url = "https://tds.hycom.org/thredds/dodsC/GLBy0.08/expt_93.0"
ds = xr.open_dataset(url, drop_variables='tau')
# ds.nbytes / 1e9 #123944.239138864
# ds

# Subset the area surrounding the transect 
# (Do I even need to do this?)
lons = [min(glon)-.1, max(glon)+.1]
lats = [min(glat)-.1, max(glat)+.1]

# only want to look temperature and salinity in 
# the upper 1000 meters.
tds = ds[['water_temp', 'salinity']].sel(
    time=slice(min(gtime), max(gtime)),
    depth=slice(0, 1000), 
    lon=slice(lons[0], lons[1]), 
    lat=slice(lats[0], lats[1])
)
# tds.nbytes / 1e9 # 1.251028936 GB
# tds

# Select the time, lon, and lat nearest the transect using advanced indexing
xds = tds.sel(
    time=xr.DataArray(gtime, dims='point'), 
    lon=xr.DataArray(glon, dims='point'),
    lat=xr.DataArray(glat, dims="point"),
    method='nearest'
)
# xds.nbytes / 1e9 # 3.432e-06
# # xds # this hangs when I try to look at the result

# I want to download a subset at each point.
# xds.to_netcdf('gofs_transect_test.nc',)