In everyday life, one almost always measures the temperature of something by bringing it to equilibrium with some calibrated system (i.e., your thermometer), and reading the temperature off of that. But as you quite rightly guess, this isn't possible with ultracold gases- they are too cold, have too few atoms, and need to be maintained in isolation from any macroscopic environment.
So instead, the temperature is determined by applying statistical mechanics. This can take various forms, depending on the temperature regime, but basically the idea is that you calculate what the momentum or spatial profile of the gas should be as a function of temperature, taking into account things like the masses of the atoms, their interactions, and the confining potential used, and then compare this prediction to actual measurements to back out the temperature.
Now, as you suspect and Mark Mitchison confirms in comments, the most common (although not only) method of measuring the density or momentum distribution of these gases is destructive, because it involves shining an on-resonant laser on the atoms that blasts them away. But this is okay for two reasons: one, because the exposure time is chosen to be short relative to the average speed of the atoms so that they mostly won't move too much, and two, because any atoms that are really speeding away are Doppler shifted off-resonance from the imaging laser and therefore don't scatter much light. So if properly set up this gives you the density distribution to a reasonable accuracy, and if you let the atoms freely expand for a time this density distribution corresponds to the in-trap momentum distribution.
As an example of how this all works, near $Tc$ there will be a significant portion of the gas that is still not in the condensate, and forms a broad thermal background. Images of the gas after an expansion period will show a density distribution like the famous Nobel pictures:
For the left two panels, which are at high and intermediate temperature, you can see the broad background which is almost completely replaced by the narrow condensate in the rightmost picture. This broad background just follows the usual Maxwell distribution of velocities, so if you can see it you can use it to determine the temperature (ignoring the central spike of the condensate). However, far below $T_C$ it will be too small to detect, and it is necessary to determine the temperature by another means.