I've figured out a way of solving it, I'm not sure if it works. Anyway here I go:
The initial temperature distribution can be written as an infinite sum of spherical functions,
T(\theta, \phi, t=0) = \sum_{l=0}^\infty \sum_{m=-l}^{+l}a_{l,m} Y_{l,m}(\theta, \phi)
The heat equation for one of those spherical functions can be solved easily because they have the property \Delta_{S^2} Y_{l,m}(\theta,\phi)=l(l+1)Y_{l,m}(\theta,\phi) with the spherical Laplace operator \Delta_{S^2} (it might be l(l-1) instead of l(l+1), I'm not sure anymore).
The solution to this particular initial condition is (the constants set equal to 1) T(\theta,\phi,t)=Y_{l,m}(\theta,\phi)e^{-l(l+1)t}.
Since the heat equation is linear one can solve it separately for each part of the infinite sum, and therefore the solution with general initial conditions is
T(\theta, \phi, t) = \sum_{l=0}^\infty e^{-l(l+1)t} \sum_{m=-l}^{+l}a_{l,m} Y_{l,m}(\theta, \phi).
So the problem is practically solved if you have found the series representation of the initial condition.
Perhaps I have saying something stupid but I think for large times the temperature would be uniform in all over the surface.
Yes, that's true.