While the zero-temperature properties of harmonically trapped cold few-atom systems have been discussed fairly extensively over the past decade, much less is known about the finite-temperature properties. Working in the canonical ensemble, we characterize small harmonically trapped atomic systems as a function of the temperature using analytical and numerical techniques. We present results for the energetics, structural properties, condensate fraction, superfluid fraction, and superfluid density. Our calculations for the two-body system underline that the condensate and superfluid fractions are distinctly different quantities. Our work demonstrates that the path-integral Monte Carlo method yields reliable results for bosonic and fermionic systems over a wide temperature range, including the regime where the de Broglie wavelength is large, i.e., where the statistics plays an important role. The regime where the Fermi sign problem leads to reasonably large signal-to-noise ratios is mapped out for selected parameter combinations. Our calculations for bosons focus on the unitary regime, where the physics is expected to be governed by the three-body parameter. If the three-body parameter is large compared to the inverse of the harmonic oscillator length, we find that the bosons form a droplet at low temperature and behave approximately like a noninteracting Bose and eventually Boltzmann gas at high temperature. The change of the behavior occurs over a fairly narrow temperature range. A simple model that reproduces the key aspects of the phase-transition-like feature, which can potentially be observed in cold atom Bose gas experiments, is presented.