Spherical harmonics popularity. I know that any periodic function of it's two arguments can be expended on spherical harmonics, and in here they claim  

The development in spherical harmonics is the equivalent, applied to the angular functions, of the development in Fourier series for the periodic functions.

I thought that this is why they are so popular in physics and mathematics but then i realized ( and correct if i am wrong ) we can simply do a Fourier expansion on a sphere.
My questions are 
What properties makes them so popular in physics and mathematics ?
Why do we usual prefer to expend in harmonics rather than Fourier ?
Is it possible to define harmonics on a circle rather than a sphere ? 
 A: For simplicity, let's talk about axisymmetric functions: functions that are only functions of the colatitude $\theta$; we can extend to more general functions once we know what's happening there.
The integration measure on a sphere for axisymmetric functions is $2\pi\sin{\theta} \, d\theta$ (having integrated out the $\phi$ part). Why? It's like polar coordinates, but the length of a circle a (spherical) distance $\theta$ from the origin at the "north pole" is $2\pi \sin{\theta}$; this gives you the "weight" you have to assign to function values at $\theta$, since it covers more area if $\theta$ is close to $\pi/2$ and we are near the "equator".
With this in mind, what should orthogonal functions on the sphere look like? Let's see what happens with two typical elements of a Fourier series of a periodic function with period $\pi$. I'm going to use the exponential form because the integrals are simpler, but with enough patience you can find the problem occurs if we use the trigonometric form.
$$ \int_0^{\pi} e^{2ni\theta}e^{-2mi\theta} \sin{\theta} \, d\theta = \frac{2}{1-4(m-n)^2} $$
(you can find this by integrating by parts a couple of times, or using $\sin{\theta}=\frac{1}{2i}(e^{i\theta}-e^{-i\theta})$) Oh dear. This integral is never zero, which means that $e^{2in\theta}$ is not an orthogonal basis: nothing's orthogonal any more! This makes all the nice formulae from Fourier analysis useless, since they're based on orthogonality.
What to do instead? Since we're on $[0,\pi]$ and it's easier to work with polynomials, we can change variables to $x=\cos{\theta}$, which is essentially the "$z$-coordinate". Then the integral of a function $f(\theta)$ becomes
$$ \int_{-1}^1 f(\arccos{x}) \sqrt{1-x^2} \, dx. $$
The axisymmetric spherical harmonics are now defined as the orthogonal polynomials for the inner product using this measure:
$$ \langle f,g\rangle = \int_{-1}^1 f(x)\overline{g(x)} \sqrt{1-x^2} \, dx, $$
where we've written our functions using $x$ as the argument instead of $\arccos{x}$; since it's a bijection it makes no difference theoretically. It turns out that the polynomials given by
$$ P_n(x) = \left(\frac{d}{dx} \right)^n (x^2-1)^n $$
are an orthogonal basis for the inner-product space of square-integrable functions on $[-1,1]$ with this inner product. These, of course, are the Legendre polynomials, up to a constant factor.

What about more general spherical harmonics? One can proceed in a similar way, but finding the orthogonal functions by the same method is quite fiddly. Instead, let's look at a different motivation, that explains the "harmonics" part.
Laplace's equation is
$$ \Delta u = 0, $$
where $\Delta = \operatorname{div}\operatorname{grad} $. I think it's reasonable to say that Laplace's equation and its relatives are quite important in physics. Let's look for homogeneous solutions: those of the form $u(r,\theta,\phi) = r^n f(\theta,\phi)$. Applying the Laplacian in spherical coordinates gives
$$ \Delta u = n(n+1) r^n f(\theta,\phi) + r^n\Delta_S f(\theta,\phi), $$
where $\Delta_S$ is the angular part of the Laplacian. Therefore $f$ satisfies
$$ \Delta_S f = -n(n+1)f, $$
and any smooth $f$ that satisfies this equation is called a spherical harmonic of weight $n$. The "angular momentum" derivative $-i\partial/\partial\phi$ commutes with $\Delta_S$, so there are simultaneous eigenfunctions of both. It turns out there are $2n+1$ of these, with "angular momenta" $-n,-n+1,\dotsc,n$. But this doesn't get us any closer to finding them, so let us tack and look at them in a different way.
We know plenty of smooth functions that satisfy Laplace's equation: there are polynomials of any degree that satisfy it: harmonic polynomials. Asking for simultaneous eigenfunctions for $\Delta$ and $-i\partial\phi$ leads to considering polynomials of the form $(x \pm iy)^{|m|} r^{n-|m|} p(z/r)$, where $p$ has degree $n-\lvert m \rvert$. Inserting this into Laplace's equation implies that $p$ must satisfy the associated Legendre equation, as for spherical harmonics derive the "usual" way. So spherical harmonics are just polynomials in $x,y,z$, chosen to be eigenfunctions of the Laplacian and $-i\partial_{\phi}$.
A: 
What properties makes them so popular in physics and mathematics ?

Your quote basically says this: they are an orthonormal basis in two variables which also turn out to be angle-space eigenfunctions for a variety of 'spherical' operators (e.g. the spherical Laplacian).

Why do we usual prefer to expend in harmonics rather than Fourier ?

It depends on the problem. Multivariable Fourier series are less frequently found as eigenfunctions of operators with physical significance, so in physics you see more spherical harmonics. In image processing you see more multidimensional Fourier series.

Is it possible to define harmonics on a circle rather than a sphere ? 

Yes. In fact it is just a Fourier series in a single variable. 
