Calculating flux integral in spherical coordinates Let $S$ be the surface of a unit sphere. 
Let $\vec{A}=\frac{1}{r^2}\hat{e}_{r}$
The flux through the surface $S$ is given by: $\int_{S}\vec{A}\space\cdot\space d\vec{S}$
$$d\vec{S}=r^2\sin\theta \,d\theta d\phi\hat{e}_{r}$$
$$\int_{S}\vec{A}\space\cdot\space d\vec{S}=\int_{s}(\frac{1}{r^2}\hat{e}_{r})\cdot(r^2\sin\theta \,d\theta d\phi\hat{e}_{r})=\int_{S}\sin\theta \space d\theta d\phi=\int_{0}^{2\pi}\int_{0}^{\pi}\sin\theta \ d\theta d\phi=4\pi$$
The above steps are all correct. 
My question is, how is this $d\vec{S}$ calculated? 
Calculating flux using Cartesian coordinates is very familiar to me, however I can't grasp how it is done in spherical coordinates or in any other general coordinate system. Especially, "the unit vector in spherical coordinates is a function of position" part is quite confusing.
I could, of course, solve the above problem by converting everything into Cartesian coordinates. But this seems like going the "long way around" when everything is already set up so nicely in spherical coordinates. 
 A: There is a handy and intuitive way to derive the surface element
$$
dS=r^2\sin\theta \,d\theta d\phi\
$$
Think of the surface area $dS$ as an infinitesimal square. In spherical coordinates, one of its sides is simply the arc length $r d\theta$ along the $\theta$-direction and the other side along the $\phi$-direction is the arc length $\rho d\phi$, where the arc lies on the circle with its radius given by $\rho=r\sin\theta$. Thus, together, one has
$$
dS= (r d\theta)( \rho d\phi)=r^2\sin\theta \,d\theta d\phi\
$$
Since the surface is on a sphere, the corresponding vector form is 
$$d\vec{S}=r^2\sin\theta \,d\theta d\phi\hat{e}_{r}$$
A: For any elemental vector area
$$
d\vec{S}=dS\,\hat{n}
$$
where $\hat{n}$ is the outward unit normal.
In spherical polars, it's a standard result (check any book discussing them) that the element of scalar area is
$$
dS=r^2\sin\theta \,d\theta d\phi\
$$
I your case $\hat{n}$ is the radial unit vector $\hat{e}_{r}$.
In addition, there is not one unit vector in spherical polar coordinates, but three, $\hat{e}_{r}, \hat{e}_{\theta}$ and $\hat{e}_{\phi}$. Three are required for the same reason you  need 3 cartesian unit vectors. You are correct that they change (direction) depending on the point in space you are considering, but remain mutually orthogonal to each other. The variation with location is part of their power when considering problems with spherical symmetry.
A: I learned that calculating directly in spherical coordinate is impractical. 
Whenever doing vector calculus, all vectors should always be converted to Cartesian, which makes the calculation magnificently simpler(due to the position independent nature of Cartesian unit vectors).
Directly coming up with a spherical expression for $d\vec{S}$ can be done on simple surfaces, such as the surface of a sphere, but for general surfaces, it can get very labor intensive. 
A: As you have stated in your comments, if you are using Cartesian coordinates, the general 'surface vector' (that is, the normal unit vector times the surface element is given by)
$$ \mathrm{d} \vec{S} = \frac{\mathrm{d} \vec{r}}{\mathrm{d} s} \times \frac{\mathrm{d} \vec{r}}{\mathrm{d} t} \mathrm{d} s \mathrm{d} t, $$
with your surface being parametrized using $s$ and $t$, that is, $S=S(s,t)$.
This equation is also true if you are working in spherical coordinates, or any curvilinear coordinate system in general. The difference is in how you take the derivative. In Cartesian coordinates, this is probably obvious, since you have, for the differential element of $\vec{r} = r_x \hat{x} + r_y \hat{y} + r_z \hat{z}$
$$ \mathrm{d} \vec{r} = \frac{\partial r_x}{\partial x}\mathrm{d} \vec{x} +  \frac{\partial r_y}{\partial y}\mathrm{d} \vec{y} +  \frac{\partial r_z}{\partial z}\mathrm{d} \vec{z}$$
However, this formula changes when you move to more general curvilinear systems. If your coordinates are still orthonormal (and spherical coordinates are orthonormal), then you will have the scale factors $h_q$ in your formula. For example in spherical coordinates
$$ \mathrm{d} \vec{r} = h_r \frac{\partial r_r}{\partial r} \mathrm{d} \vec{r} +   h_{\theta} \frac{\partial r_{\theta}}{\partial \theta} \mathrm{d} \vec{\theta} + h_{\phi} \frac{\partial r_{\phi}}{\partial \phi} \mathrm{d} \vec{\phi}$$
For example you could ask what the velocity is in spherical coordinates, and you would find that the general expression is
$$ \frac{\mathrm{d} \vec{r}}{\mathrm{d} t} = \dot{r} \hat {r} + r \sin (\phi)  \dot{\theta} \hat{\theta} +  r \dot{\phi} \hat{\phi} $$
Now, since it doesn't matter if you take the derivative with respect to time $t$, or an arbitrary parameter $t$ or $s$, the above formula gives you the expression on how to correctly take derivatives in spherical coordinates. If you do this consistently with your parametrization, then evaluate the cross product with this result, then your surface element will be properly scaled.
