Question: I know that both Fourier and wavelet can be used for compression of signals. 

The Fourier series guarantees that it gets the closest approximation of the original signal in the least squared sense. It also decomposes the signal into orthogonal components. So, I wonder why it isn't obvious to always use the `Fourier` series as a compression algorithm?

The Wavelet transform allows us to observe frequencies at different times, because it is localized in time (at least relative to the fourier transform). So it is only useful when we care about the interpretability of the components, it seems to me.

Why is the `wavelet compression` used in practice much more than compression using a `fourier series`?

Answer: 
In short: **wavelets and relatives are pretty good at compacting a sufficiently large class of regular-enough and useful signals and images**. What follows are signal properties (and corresponding wavelet features) that make wavelets good (not best) candidates for lossy compression:

* **piecewise-smooth** (*vanishing moments*, or cancellation of low-order polynomials)
* **edges or jumps** (*gradient*-like or Laplacian *behavior* of wavelets)
* **localized oscillations** (zero-average and *wiggling* wavelet *shape*)
* **noise or spurious events** (*orthogonality* and sparsity enhancement)

As said by @hops, the efficiency of wavelets for compression depends of the good matching of a signal class and the chosen wavelets.

Let us restrict here to non-redundant discrete transformations: discrete Fourier and discrete wavelets. Both are orthogonal, or close enough (bio-orthogonal wavelets) to skip the distinction. 

So both, when transform coefficients are discarded, are least-squares approximations. But least-squares are not the best metric for compression: if you double the amplitude of a sample, the energy is multiplied by $4$, but it only adds one bit to the stored data, in a $\log\_2$ reasoning.

In a way, a transform will help compression if it reduces the logarithmic cost or bit-budget of a discrete signal; hence, if the coefficients have a power law $1/c^{\beta}$ that decreases fast enough (the highest the $\beta$, the better).

This feature is often called the **compressibility of signals**. It can be assessed empirically, but also theoretically based on complicated functional analysis (Besov, Sobolev spaces). 

However, consider the useful class of $C^\alpha$ piecewise regular signals, $\alpha \ge 1$. They are locally smooth, with jumps (edges). Taking the largest $M$ coefficients for compression (or **nonlinear** **approximation**), the mean squared error for the Fourier basis will asymptotically decrease in $1/M$ for Fourier (because of Gibbs ringing). While wavelet approximations can decrease as $1/M^{-2\alpha}$. **The smoother the signal, the better the wavelet compressibility**. In other words, Fourier cannot use the regularity of the signal.

In practice, this is not so simple. The complex aspect of Fourier is somewhat difficult to handle. Quantization should be taken into account. The storage of the highest coefficients location costs too. Perceptual distortion should play a role, at least for lossy compression. 

So at low compression rates, local Fourier decompositions like the JPEG DCT can perform as well as wavelets, since asymptotic proofs do not apply. Indeed, local cosine bases, lapped orthogonal transforms that bridge the gap between Fourier and standard wavelets perform well too, and they are used in the MP3 standard (MDCT). For images, the 2D aspect renders separable wavelets less efficient than in 1D.

For strict lossless compression through, (wavelet) transforms are not the state-of-the-art yet.
Fourier analysis is great for signals with unchanging characteristics. So it is usually a good fit for modeling noise. Unfortunately, most of the time we are not interested in noise but in signals. You can try regaining a signal by subtracting noise, but straightforward subtraction requires the *phase* of the noise (which does not lead itself well to modeling for stochastic signals), and spectral subtraction is basically selective attenuation that kills weak parts of the signal along with the noise.

So one usually wants to model the signal more than the noise, and the whole point of a "signal" is its information content, namely intentionally *changing* signal characteristics. For that, local support of your basis functions is a reasonable approach since it helps keeping the nonzero coefficients for representing the changing signal better localized in space-time.