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}$.
