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.
Best Answer
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.