Most analysis of Cosmic Microwave Background spherical harmonic coefficients a_lm has focused on estimating the power spectrum C_l=<|a_lm|^2> rather than the coefficients themselves. We present a minimum-variance method for measuring a_lm given anisotropic noise, incomplete sky coverage and foreground contamination, and apply it to the WMAP data. Our method is shown to constitute lossless data compression in the sense that the widely used quadratic estimators of the power spectrum C_l can be computed directly from our a_lm-estimators. As the Galactic cut is increased, the error bars Delta-a_lm on low multipoles go from being dominated by foregrounds to being dominated by sample variance from other multipoles, with the intervening minimum defining the optimal cut. Applying our method to the WMAP quadrupole and octopole, we find that their previously reported "axis of evil" alignment appears to be rather robust to Galactic cut and foreground contamination.