When you do an integral over a surface in coordinates, we use basis vectors to talk about directions and their cross products to talk about areas. For instance, we can parameterize a surface in (edit) $s, t$ and talk about the differential area $d\mathbf A= ( \partial_s \mathbf r \times \partial_t \mathbf r) \, ds \, dt$. Of course, if the surface if one of constant $z$, then we can parameterize in terms of $x,y$ and $\partial_x \mathbf r = \hat{\mathbf x}$, and $\partial_y \mathbf r = \hat{\mathbf y}$, so we get $d\mathbf A = \hat{\mathbf x} \times \hat{\mathbf y} \, dx \, dy$.
In (15), the same approach is used, but using $\theta, \phi$ instead of $x, y$. You could use the Jacobian and get the same answer by expressing $dx$ in terms of $dr, d\theta, d\phi$ and similarly for $dy$, but it's not really necessary.
In (23), it's not clear to me why this integral is phrased in terms of $r$ at all, since the surface is a sphere and thus lies at constant $r = 2$. Nevertheless, you can use the same approach.
$$d\mathbf A = (\partial_r \mathbf r \times \partial_\theta \mathbf r) \, dr \, d\theta$$
We get $\partial_r \mathbf r = \hat{\mathbf r}$ and $\partial_\theta \mathbf r = r \hat{\boldsymbol \theta}$. That explains where the extra $r$ comes from.
(47) is worked the same way, again using $\theta, \phi$.
Now, you may be wondering, if you can find all these area elements without using the Jacobian, then why do people say "use the Jacobian"? Well, that's because, if you already have the vector of the area element expressed in one coordinate system, you can move it to another using the Jacobian. Here, we don't have the area element vector expressed in a coordinate system yet, so it doesn't make sense to use (say) Cartesian and then push it forward with the Jacobian.
It's easiest to demonstrate this with a line integral. Consider a line integral with differential $d\boldsymbol \ell = \partial_x \mathbf r \, dx$. Let $u = u(x)$ be some function, so that the integral can be equivalently phrased as $d\boldsymbol \ell = (\partial_u x)(\partial_x \mathbf r) \, du = \partial_u \mathbf r \, du$. This all comes from the chain rule, and $\partial_u x$ is the 1d version of the Jacobian.
So what are we doing with the Jacobian? We're taking derivatives of $\mathbf r$ with respect to one set of coordinate and converting them to derivatives with respect to another set of coordinates. The solutions here skip the Jacobian because they take derivatives in spherical coordinates directly.
Edit:
1) You're correct; I was talking about a surface of constant $z$ (and was lazy in doing so).
2) My point is merely that you can derive $d\mathbf A = r \hat{\boldsymbol \phi} \, dr \, d\theta$ quite easily in a direct fashion, rather than expressing the vector in cartesian, writing out the Jacobian matrix, and then pushing it forward into a spherical coordinate system.
Since they know the answer, I think honestly they just knew $dA = r \, dr \, d\theta$ and thought nothing of it. After working with spherical and polar coordinate systems long enough, you kind of just know these things. But I must emphasize that that, in itself, is not using the Jacobian matrix, because the Jacobian matrix must act on some vector, converting it from one coordinate system to another. Since they don't write down the cartesian components of $\hat{\mathbf n}$, it doesn't seem to me that they're actually using it.
3) I'm referring to something more fundamental. They say quite clearly that this should be "the surface of a sphere", which by definition doesn't vary with $r$, and so $r$ should not be integrated over, and yet they end up with an expression that does depend on $r$, which makes no sense. It's like they're integrating over some completely different surface from what they describe.
As for using spherical coordinates, yeah, if you're just going to transform into those coordinates anyway, they should be used from the outset. Parameterizing these surfaces in another coordinate system is largely a pointless exercise; one of the main advantages of working in other coordinate systems is to simplify integrals and to match the natural symmetries of a given problem.
Best Answer
hint
Use the formula
$$\iint_S\vec{F}(x,y,z).\vec{ds}=\iiint_Vdiv(\vec{F})dv$$
to get $$2\iiint_V(x+y+z)dxdydz=$$ $$2\int_0^1\int_0^{2\pi}\int_0^{\frac{\pi}{2}}(\cos(t)\sin(f)+\sin(t)\sin(f)+\cos(f))\sin(f)drdtdf$$