Surface integrals in spherical coordinates If I am given a surface in spherical coordinates $(r,\theta,\varphi)$, such that it is parametrised as:
$$
\begin{align}
r&=r(\theta,\varphi)\\
\theta&=\theta\\
\varphi&=\varphi
\end{align}
$$
What is the area $S$ of such surface? Or more specifically, can you show how to get the result:
$$
S=\int_{0}^{2\pi}\int_{0}^{\pi}\sqrt{r^2+\left(\frac{\partial r}{\partial \theta}\right)^2 + \frac{1}{\sin^2\theta}\left(\frac{\partial r}{\partial \varphi}\right)^2}\;r\sin\theta\;{\rm d}\theta\,{\rm d}\varphi
$$

Some definitions that I am using:
$k$-surface:
Let $k,N\in\mathbb{N}$, $k<N$, $M\subset \mathbb{R}^N$ is called a $k$-surface, if there exists a non-empty open set $E\subset \mathbb{R}^k$ and a map $\varphi:\mathbb{R}^k\to \mathbb{R}^N$, such that: (i) $\varphi(E)=M$, (ii) $\varphi\in C^1(E;\mathbb{R}^N)$, and (iii) the rank of Jacobi matrix of $\varphi$ is equal $k$ everywhere on $E$.
The surface is called simple if $\varphi$ is also injective on $E$ and $\varphi^{-1}$ is continuous of $\varphi(E)$.
Surface integral of the first kind:
Let $k,N\in\mathbb{N}$, $k<N$, $M\subset \mathbb{R}^N$ is a simple $k$-surface parametrized by the map $\varphi:\mathbb{R}^k\to \mathbb{R}^N$, from the open set $E\subset \mathbb{R}^k$ and $f:\mathbb{R}^N\to\mathbb{R}$ is defined on $M$. The surface integral of the first kind is defined by:
$$
\int_M f\,\mathrm{d}S:=\int_E f(\varphi(t))\sqrt{\det{G(D_\varphi(t))}}\,\mathrm{d}t\,,
$$
if the integral on the right exists in the Lebesgue sense and is finite. Here, $G(A)$ denotes the Gramm matrix made from columns of $A$ and $D_\varphi$ is the Jacobi matrix of the map $\varphi$. The numeric value of:
$$
S_k(M):=\int_M f\,\mathrm{d}S\,,
$$
is called the $k$-dimensional surface area of the $k$-surface $M$.
Motivation for the question
Now these definitions can be used to calculate e.g. the surface of a unit sphere. When one describes the sphere by the map (omitting one longitudinal line): $\varphi: (\eta,\psi)\mapsto(\cos\psi\cos\eta,\cos\psi\sin\eta,\sin\psi)$, where $(\eta,\psi)\in E=(0,2\pi)\times(-\frac{\pi}{2},\frac{\pi}{2})$. The Gramm matrix looks like:
$$
\begin{pmatrix}
\cos^2\psi & 0 \\
0 & 1
\end{pmatrix}
$$
one ends up with the (here $\lambda_2$ denotes the Lebesgue measure):
$$
S_2(M) = \int_E 1\sqrt{\det{G(D_\varphi(\eta,\psi))}}\,\mathrm{d}\lambda_2(\eta,\psi)=\int_0^{2\pi}\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos\psi\,\mathrm{d}\psi\,\mathrm{d}{\eta}=4\pi
$$
I did not do anything else but blindly followed the definitions.
The wrong result:
If I do the same approach to the above problem at hand. My map is $\varphi:\,(\theta,\phi)\mapsto (r(\theta,\phi),\theta,\phi)$, the Jacobian is then:
$$
\begin{pmatrix}
\frac{\partial r(\theta,\phi)}{\partial \theta} & \frac{\partial r(\theta,\phi)}{\partial \phi} \\
1 & 0 \\
0 & 1
\end{pmatrix}
$$
The rank of this matrix is 2 as needed for a 2-surface. The Gramm matrix is then:
$$G(D_\varphi)=
\begin{pmatrix}
1+ \left(\frac{\partial r(\theta,\phi)}{\partial \theta}\right)^2 & \frac{\partial r(\theta,\phi)}{\partial \theta}\frac{\partial r(\theta,\phi)}{\partial \phi} \\
\frac{\partial r(\theta,\phi)}{\partial \theta}\frac{\partial r(\theta,\phi)}{\partial \phi} & 1+ \left(\frac{\partial r(\theta,\phi)}{\partial \phi}\right)^2 
\end{pmatrix}
$$
with determinant evaluated to:
$$
\det(G(D_\varphi)) = 1 + \left(\frac{\partial r(\theta,\phi)}{\partial \theta}\right)^2 + \left(\frac{\partial r(\theta,\phi)}{\partial \phi}\right)^2
$$
which leads to the surface area:
$$S=\int_{S\subset\mathcal{R}^3}{\rm d}S=\int_{0}^{2\pi}\int_{0}^{\pi}\sqrt{1+\left(\frac{\partial r}{\partial \theta}\right)^2 + \left(\frac{\partial r}{\partial \varphi}\right)^2}\;{\rm d}\theta\,{\rm d}\varphi$$
Which is incorrect. The question is why, and please show the correct way with explanation because as suggested by the example with unit sphere, I did not have to do any transformations and it did work with curvilinear coordinates right away. So the answer provided by Quanto does not address this at all.
Note to the edit:
I have added quite a lot of details since the answer by Quanto, but the root problem is the same - to see by calculation explicitly how one arrives to the correct result and to understand why the my calculation for a sphere works (where i am not using Cartesian coordinates) but fails here.
 A: The map $\varphi:(\theta,\phi)\rightarrow\Big(f(\theta,\phi),\theta,\phi\Big)$ is not a parametric representation of the spherical surface $r=f(\theta,\phi)$. Rather, $\varphi$ is a parametric representation of the surface $x=f(y,z)$ since parametric surfaces (and curves) are inherently represented in Cartesian form. What you're doing is equivalent to saying $\theta \longrightarrow \Big(\theta^2,\theta\Big)$ is a parametric represention of the polar spiral $r=\theta^2$ when, in reality, $\theta \rightarrow \Big(\theta^2,\theta\Big)$ is parametric representation of the parabola $x=y^2$ while $\theta \rightarrow \Big(\theta^2 \cos(\theta),\theta^2 \sin(\theta)\Big)$ is the spiral $r=\theta^2$.
Let's refer to your example in which you seek to compute the surface area of a sphere.
In spherical coordinates, the equation of a sphere is $r=1$ on the domain $(\theta,\phi)\in[0,2\pi)\times [0,\pi]$. You can represent this parametrically as $$(\phi,\theta) \longrightarrow \Big(\sin(\phi)\cos(\theta),\sin(\phi)\sin(\theta),\cos(\phi)\Big)$$ simply by converting from spherical to cartesian coordinates. However, the image of your function $\varphi:(\phi,\theta)\longrightarrow (1,\phi,\theta)$ on the same domain is the rectangular patch $\{1\}\times [0,2\pi)\times [0,\pi]$ embedded in the vertical plane $x=1$ which has area $2\pi^2$. This is precisely why $$\int_0^{2\pi}\int_0^{\pi}\sqrt{1+\big(f_{\theta}\big)^2+\big(f_{\phi}\big)^2}d\phi d\theta=2\pi^2 \neq 4\pi$$
whenever $f(\phi,\theta)=1$; this integral is calculating the area of the surface $x=1$ on $(y,z)\in [0,2\pi)\times[0,\pi]$ which is the image of your map $\varphi$.
If you want to find the area of the surface given in spherical coordinates by $r=f(\phi,\theta)$ defined on domain  $(\phi,\theta)\in \mathcal{U}$ you will need to express this parametrically as $$\vec{p}(\phi,\theta)=\Big(f(\phi,\theta)\sin(\phi)\cos(\theta),f(\phi,\theta)\sin(\phi)\sin(\theta),f(\phi,\theta)\cos(\phi)\Big)$$ The area of the surface will be $$S=\int \int _{\mathcal{U}}||\vec{p}_\phi \times \vec{p}_\theta||d\phi d\theta$$ If you compute $||\vec{p}_\phi \times \vec{p}_\theta||$ you will surely obtain your desired result.
A: The standard surface integral, or the first approach, is for the cartesian coordinates, i.e.,
$$S=\int_{S\subset\mathcal{R}^3}{\rm d}S=\int\int \sqrt{1+\left(\frac{\partial z}{\partial x}\right)^2 + \left(\frac{\partial z}{\partial y}\right)^2}\;{\rm d}x\,{\rm d}y$$
It can not be simply recast into one with spherical coordinates.
A: There is an elementary way to get the areal element, but it's kind of tedious. Along the surface let
$$\vec r=\langle x,y,z\rangle=\langle r\sin\theta\cos\phi,r\sin\theta\sin\phi,r\cos\theta\rangle$$
Then
$$\begin{align}d\vec r&=\left\langle\frac{\partial r}{\partial\theta}\sin\theta\cos\phi+r\cos\theta\cos\phi,\frac{\partial r}{\partial\theta}\sin\theta\sin\phi+r\cos\theta\sin\phi,\frac{\partial r}{\partial\theta}\cos\theta-r\sin\theta\right\rangle\,d\theta\\
&\quad+\left\langle\frac{\partial r}{\partial\phi}\sin\theta\cos\phi-r\sin\theta\sin\phi,\frac{\partial r}{\partial\phi}\sin\theta\sin\phi+r\sin\theta\cos\phi,\frac{\partial r}{\partial\phi}\cos\theta\right\rangle\,d\phi\end{align}$$
This writes an infinitesimal displacement along the surface in terms of changes in the independent variables. We find a vector normal to the surface whose magnitude is the area of a parallelogram two of whose sides are the two infinitesimal vectors above:
$$\begin{align}d^2\vec A&=\pm\left\langle\frac{\partial r}{\partial\theta}\sin\theta\cos\phi+r\cos\theta\cos\phi,\frac{\partial r}{\partial\theta}\sin\theta\sin\phi+r\cos\theta\sin\phi,\frac{\partial r}{\partial\theta}\cos\theta-r\sin\theta\right\rangle\,d\theta\\
&\quad\times\left\langle\frac{\partial r}{\partial\phi}\sin\theta\cos\phi-r\sin\theta\sin\phi,\frac{\partial r}{\partial\phi}\sin\theta\sin\phi+r\sin\theta\cos\phi,\frac{\partial r}{\partial\phi}\cos\theta\right\rangle\,d\phi\\
&=\pm\left\langle\frac{\partial r}{\partial\theta}\frac{\partial r}{\partial\phi}\sin\theta\cos\theta\sin\phi+r\frac{\partial r}{\partial\phi}\cos^2\theta\sin\phi\right.\\
&\quad-\frac{\partial r}{\partial\theta}\frac{\partial r}{\partial\phi}\sin\theta\cos\theta\sin\phi+r\frac{\partial r}{\partial\phi}\sin^2\theta\sin\phi-r\frac{\partial r}{\partial\theta}\sin\theta\cos\theta\cos\phi+r^2\sin^2\theta\cos\phi,\\
&\quad\frac{\partial r}{\partial\theta}\frac{\partial r}{\partial\phi}\sin\theta\cos\theta\cos\phi-r\frac{\partial r}{\partial\phi}\sin^2\theta\cos\phi-r\frac{\partial r}{\partial\theta}\sin\theta\cos\theta\sin\phi+r^2\sin^2\theta\sin\phi\\
&\quad-\frac{\partial r}{\partial\theta}\frac{\partial r}{\partial\phi}\sin\theta\cos\theta\cos\phi-r\frac{\partial r}{\partial\phi}\cos^2\theta\cos\phi,\\
&\quad\frac{\partial r}{\partial\theta}\frac{\partial r}{\partial\phi}\sin^2\theta\sin\phi\cos\phi+r\frac{\partial r}{\partial\phi}\sin\theta\cos\theta\sin\phi\cos\phi+r\frac{\partial r}{\partial\theta}\sin^2\theta\cos^2\phi\\
&\quad+r^2\sin\theta\cos\theta\cos^2\phi-\frac{\partial r}{\partial\theta}\frac{\partial r}{\partial\phi}\sin^2\theta\sin\phi\cos\phi-r\frac{\partial r}{\partial\phi}\sin\theta\cos\theta\sin\phi\cos\phi\\
&\quad\left.+r\frac{\partial r}{\partial\theta}\sin^2\theta\sin^2\phi+r^2\sin\theta\cos\theta\sin^2\phi\right\rangle\,d\theta\,d\phi\\
&=\pm\left\langle r\frac{\partial r}{\partial\phi}\sin\phi-r\frac{\partial r}{\partial\theta}\sin\theta\cos\theta\cos\phi+r^2\sin^2\theta\cos\phi,\right.\\
&\quad-r\frac{\partial r}{\partial\phi}\cos\phi-r\frac{\partial r}{\partial\theta}\sin\theta\cos\theta\sin\phi+r^2\sin^2\theta\sin\phi,\\
&\quad\left.r\frac{\partial r}{\partial\theta}\sin^2\theta+r^2\sin\theta\cos\theta\right\rangle\,d\theta\,d\phi\end{align}$$
So now
$$\begin{align}d^2A&=\left\lVert d^2\vec A\right\rVert=\left\{r^2\left(\frac{\partial r}{\partial\phi}\right)^2\sin^2\phi+r^2\left(\frac{\partial r}{\partial\theta}\right)^2\sin^2\theta\cos^2\theta\cos^2\phi+r^4\sin^4\theta\cos^2\phi\right.\\
&\quad+r^2\left(\frac{\partial r}{\partial\phi}\right)^2\cos^2\phi+r^2\left(\frac{\partial r}{\partial\theta}\right)^2\sin^2\theta\cos^2\theta\sin^2\phi+r^4\sin^4\theta\sin^2\phi\\
&\quad-2r^2\left(\frac{\partial r}{\partial\phi}\right)\left(\frac{\partial r}{\partial\theta}\right)\sin\theta\cos\theta\sin\phi\cos\phi+2r^3\left(\frac{\partial r}{\partial\phi}\right)\sin^2\theta\sin\phi\cos\phi\\
&\quad-2r^3\left(\frac{\partial r}{\partial\theta}\right)\sin^3\theta\cos\theta\cos^2\phi+2r^2\left(\frac{\partial r}{\partial\phi}\right)\left(\frac{\partial r}{\partial\theta}\right)\sin\theta\cos\theta\sin\phi\cos\phi\\
&\quad-2r^3\left(\frac{\partial r}{\partial\phi}\right)\sin^2\theta\sin\phi\cos\phi-2r^3\left(\frac{\partial r}{\partial\theta}\right)\sin^3\theta\cos\theta\sin^2\phi\\
&\quad\left.+r^2\left(\frac{\partial r}{\partial\theta}\right)^2\sin^4\theta+2r^3\left(\frac{\partial r}{\partial\theta}\right)\sin^3\theta\cos\theta+r^4\sin^2\theta\cos^2\theta\right\}^{1/2}\,d\theta\,d\phi\\
&=\left\{r^2\left(\frac{\partial r}{\partial\phi}\right)^2+r^2\left(\frac{\partial r}{\partial\theta}\right)^2\sin^2\theta\cos^2\theta+r^4\sin^4\theta-2r^3\left(\frac{\partial r}{\partial\theta}\right)\sin^3\theta\cos\theta\right.\\
&\quad\left.+r^2\left(\frac{\partial r}{\partial\theta}\right)^2\sin^4\theta+2r^3\left(\frac{\partial r}{\partial\theta}\right)\sin^3\theta\cos\theta+r^4\sin^2\theta\cos^2\theta\right\}^{1/2}\,d\theta\,d\phi\\
&=\left\{r^2\left(\frac{\partial r}{\partial\phi}\right)^2+r^2\left(\frac{\partial r}{\partial\theta}\right)^2\sin^2\theta+r^4\sin^2\theta\right\}^{1/2}\,d\theta\,d\phi\\
&=\sqrt{\frac1{\sin^2\theta}\left(\frac{\partial r}{\partial\phi}\right)^2+\left(\frac{\partial r}{\partial\theta}\right)^2+r^2}\,r\sin\theta\,d\theta\,d\phi\end{align}$$
A lot of typing but we have the correct answer if the typo will. The areal element given in the original question was wrong at least because quantities with different units ($\text{length}^2$ and no units) were being added under the square root.
EDIT: There is a much more efficient way to organize the arithmetic. Since
$$\vec r=\langle r\sin\theta\cos\phi,r\sin\theta\sin\phi,r\cos\theta\rangle$$
we have
$$\begin{align}d\vec r&=\langle\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta\rangle\,dr+\langle r\cos\theta\cos\phi,r\cos\theta\sin\phi,-r\sin\theta\rangle\,d\theta\\
&\quad+\langle-r\sin\theta\sin\phi,r\sin\theta\cos\phi,0\rangle\,d\phi\\
&=\langle\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta\rangle\,dr+\langle \cos\theta\cos\phi,\cos\theta\sin\phi,-\sin\theta\rangle\,rd\theta\\
&\quad+\langle-\sin\phi,\cos\phi,0\rangle\,r\sin\theta d\phi\\
&=\hat rdr+\hat{\theta}rd\theta+\hat{\phi}r\sin\theta d\phi\end{align}$$
Since we are considering $\theta$ and $\phi$ to be the independent variables,
$$dr=\left(\frac{\partial r}{\partial\theta}\right)d\theta+\left(\frac{\partial r}{\partial\phi}\right)d\phi=\frac1r\left(\frac{\partial r}{\partial\theta}\right)rd\theta+\frac1{r\sin\theta}\left(\frac{\partial r}{\partial\phi}\right)r\sin\theta d\phi$$
So along the surface
$$d\vec r=\left(\frac1r\left(\frac{\partial r}{\partial\theta}\right)\hat r+\hat{\theta}\right)rd\theta+\left(\frac1{r\sin\theta}\left(\frac{\partial r}{\partial\phi}\right)\hat r+\hat{\phi}\right)r\sin\theta d\phi$$
Then since we can show that $\hat r\times\hat{\theta}=\hat{\phi}$, $\hat{\theta}\times\hat{\phi}=\hat r$, and $\hat{\phi}\times\hat r=\hat{\theta}$,
$$\begin{align}d^2\vec A&=\pm\left(\frac1r\left(\frac{\partial r}{\partial\theta}\right)\hat r+\hat{\theta}\right)rd\theta\times\left(\frac1{r\sin\theta}\left(\frac{\partial r}{\partial\phi}\right)\hat r+\hat{\phi}\right)r\sin\theta d\phi\\
&=\pm\left(-\frac1r\frac{\partial r}{\partial\theta}\hat{\theta}-\frac1{r\sin\theta}\frac{\partial r}{\partial\phi}\hat{\phi}+\hat r\right)r^2\sin\theta d\theta d\phi\end{align}$$
Since we can also show that $\hat r\cdot\hat r=\hat{\theta}\cdot\hat{\theta}=\hat{\phi}\cdot\hat{\phi}=1$ and that $\hat r\cdot\hat{\theta}=\hat{\theta}\cdot\hat{\phi}=\hat{\phi}\cdot\hat r=0$ we have
$$\begin{align}d^2A&=\left\lVert d^2\vec A\right\lVert\\
&=\sqrt{\left(-\frac1r\frac{\partial r}{\partial\theta}\hat{\theta}-\frac1{r\sin\theta}\frac{\partial r}{\partial\phi}\hat{\phi}+\hat r\right)\cdot\left(-\frac1r\frac{\partial r}{\partial\theta}\hat{\theta}-\frac1{r\sin\theta}\frac{\partial r}{\partial\phi}\hat{\phi}+\hat r\right)}\,r^2\sin\theta d\theta d\phi\\
&=\sqrt{\frac1{r^2}\left(\frac{\partial r}{\partial\theta}\right)^2+\frac1{r^2\sin^2\theta}\left(\frac{\partial r}{\partial\phi}\right)^2+1}\,r^2\sin\theta d\theta d\phi\end{align}$$
Much nicer :)
A: Your surface corresponds to the following idea: $$S:\quad r(\theta,\phi):=R + h(\theta,\phi)\ ,\tag{1}$$
where $R$ is the earth radius at sea level, and $(\theta,\phi)$ is the height over sea level of the point with geographical coordinates $(\theta,\phi)$. You want to know the area of $S$. In order to find this area we have to resort to cartesian coordinates, because $S$ inherits its area from the euclidean ${\mathbb R}^3$, and not from the "practical" representation. This means that we have to replace $(1)$ by
$$S:\quad (\theta,\phi)\mapsto{\bf r}(\theta,\phi)=r(\theta,\phi)\,\bigl(\cos\theta\cos\phi,\cos\theta\sin\phi,\sin\theta\bigr)\ ,\tag{2}$$
where $(\theta,\phi)\in B:=\bigl[-{\pi\over2},{\pi\over2}\bigr]\times[-\pi,\pi]$. Note that $$(x,y,z)=\bigl(\cos\theta\cos\phi,\cos\theta\sin\phi,\sin\theta\bigr)$$ is the point on the unit sphere $S^2$ having geographical coordinates $(\theta,\phi)$.
In calculus 102 it is taught that the area of the surface $(2)$ is computed as
$${\rm area}(S)=\int_B \bigl|{\bf r}_\theta\times {\bf r}_\phi\bigr|\>{\rm d}(\theta,\phi)=\int_{-\pi/2}^{\pi/2}\int_{-\pi}^\pi \Psi(\theta,\phi)\>d\phi\>d\theta\ ,$$
where $\Psi(\theta,\phi):=\bigl|{\bf r}_\theta\times {\bf r}_\phi\bigr|_{(\theta,\phi)}$
has to be computed carefully from $(2)$.
In theory one arrives at this expression $\bigl|{\bf r}_\theta\times {\bf r}_\phi\bigr|_{(\theta,\phi)}$ by noting that a tiny $[\theta,\theta+\Delta\theta]\times[\phi,\phi+\Delta\phi]$ rectangle  in the parameter plane is mapped by $(2)$ onto a tiny parallelogram of area $\bigl|{\bf r}_\theta\times {\bf r}_\phi\bigr|_{(\theta,\phi)}\>\Delta\theta\Delta\phi$.
