Here's the solution.

idx_x = xr.DataArray(myIndices[:,0], dims="points")
idx_y = xr.DataArray(myIndices[:,1], dims="points")

myValues = da.isel(x=idx_x, y=idx_y)
myValues = myValues.values # this converts it from an xarray into just an array

I went from ~35 second run times to ~0.0035 second run times. Four orders of magnitude of improvement! Woot woot!

I don't understand why putting the indices into an xarray.DataArray with dimension name "points" (rather than just as a list as I tried before) causes the .isel to work correctly, but it does.