I'm searching for an intuitive explanation of solid angles as a natural 3-dimensional analogue of angles.
It's not sound yet, but I would like to say that the length of the arc occupied by the projection of a bounded 1-dimensional hypersurface $M_1$ (i.e. a finite line segment) onto the unit circle (the $1$-sphere $S_1$) is precisely the angle in radian measure spanned by the two vectors from the origin to the "endpoints" of $M_1$ (if $M_1$ is line segment, then these endpoints are precisely the two vectors constituting the boundary $\partial M_1$ of $M_1$).
Analogously, the area occupied by the projection of a bounded 2-dimensional hypersurface $M_2$ on the unit 2-sphere $S_2$ is precisely the solid angle spanned by the "contour" of $M_2$ (which seemes, intuitively, again to be given by the boundary $\partial M_2$ of $M_2$).
As I wrote before, it's not sound yet. If someone could help me to phrase this rigorously, I'd really appreciate it.
In $2$ dimensions (for a hypersurface $M$, with one fewer for $\partial M$), a circle centred on the origin allows us to simplify $ds^2=dr^2+r^2d\theta^2$ to $ds=rd\theta$, since $dr=0$, so the arc length is $r\int d\theta$ (all integrals herein are definite, with an invisible-for-brevity integration domain). This provides a natural measure of angles in radians, first for such a circle by dividing an arc length by the radius, and then for more general arcs by defining for them a radius of curvature. We want to repeat this in $3$ dimensions. But it's not enough to derive an arc length from $ds^2=dr^2+r^2d\theta^2+r^2\sin^2\theta d\phi^2$, since the solid angle we measure in steradians requires an area $\propto r^2$ instead of an arc length $\propto r$.
So let's back up a bit. We can instead get $rd\theta$ as the arc length by dividing $dxdy$ by $dr$, to get an area accrued per increase in radius if we allow $r$ to vary. Similarly, we get $\frac{dxdydz}{dr}=r^2\sin\theta d\theta d\phi$ one dimension up. So the number of steradians, over a patch of a $2$-dimensional surface embedded in $3$-dimensional Euclidean space, is a definite integral $\int\sin\theta d\theta d\phi$. You can easily use these approaches to show a circle comprises $2\pi$ radians, while a sphere comprises $4\pi$ steradians.
There's still a subtlety, though, which arises whether we take my approach as above, or work in terms of the vectors you propose. A $1$-dimensional path in $2$-dimensional Euclidean space may bend often enough that the values of $\theta$ at its endpoints don't give its full length viz. $s=\int rd\theta$, because we need to add up the $\theta$-increasing and $\theta$-decreasing parts and give each a non-negative contribution. Writing down something like $\int r|d\theta|$ is the wrong approach! What we do is endow the path with an affine parameter increasing along it, say $u$, making the arclength $\int r|d\theta/du|du$. This may seem unimportant if your interest is only in a projection's arclength & subtended angle, but it's very important when we move up one dimension because the mathematics that then makes sense of that subtlety sheds light on what different choices of coordinates do (or don't do).
So let's apply this with one more dimension. We need something like $\int|\sin\theta||\color{blue}{\theta_u\phi_v-\theta_v\phi_u}|dudv$, where $\theta_u:=\partial\theta/\partial u$ etc. with parameters $u,\,v$ appear in the blue Jacobian determinant. This may look very different from what Wikipedia does, but it's not; in their notation (which considers an orientable surface for which we can drop some modulus signs), $\hat{r}\cdot\hat{n}=\sin\theta$ while $(\theta_u\phi_v-\theta_v\phi_u)dudv=dS/r^2$ in terms of an infinitesimal area $dS$, so both sides are dimensionless. The end result of all this is we can be confident $\int\frac{\hat{r}\cdot\hat{n}}{r^2}dS$ has the same meaning regardless of our convenience-driven choice of coordinate system.
For all this to work, we need an orientable surface (which, to be fair, a projection will be). The basic issue here is that, whereas a $1$-dimensional path need only meet quite mild conditions for us to define a direction on it and hence an affine parameter, a $2$-dimensional surface patch can be traversed in much more varied ways, and we can't use this to order the points.