As a preamble, if we have a surface given by $z=z(x,y)$ in $mathbb{R}^3$, i.e. $z =z(x,y) = sqrt{a^2-x^2-y^2}$ then $z- z(x,y)=0$ gives us a constant surface of a scalar field $f$, the $grad$ of which will tell us the normal vector to our surface: $nabla f = mathbf{k} -frac{∂z}{∂x}mathbf{i}-frac{∂z}{∂y}mathbf{j}$, the norm of which is $sqrt{1+left(frac{partial z}{partial x}right)^2+left(frac{partial z}{partial y}right)^2}$
A hemisphere of radius $a$ (with the $z$ axis as the axis of symmetry) is described by the equation $x^2+y^2+z^2 (=r^2) = a^2$ and has a surface area given by
$$iint dS = iint sqrt{1+left(frac{partial z}{partial x}right)^2+left(frac{partial z}{partial y}right)^2}dA$$
If you project onto the $xy$ plane.
This formula is obtained by considering an element of surface area $mathbf{dS}$ and its projected area $mathbf{dA}$, and relating the two by $dA = cos(alpha) dS = mathbf{hat{n}}cdotmathbf{k}, dS rightarrow dS = dA/mathbf{hat{n}}cdotmathbf{k} = sqrt{1+left(frac{partial z}{partial x}right)^2+left(frac{partial z}{partial y}right)^2}dA $ using $mathbf{n}$, the normal to the surface, from the preamble and $mathbf{k}$ the unit vector in the $z$ direction.

For the above example we can use $z =sqrt{a^2-x^2-y^2}$ to find our partial derivatives, and then evaluate the integral using plane polar co-ordinates ($0 leq theta leq 2pi, 0 leq r leq a $) to reach an answer of $2pi a^2$ as expected.
However, if we use a normal vector to this surface expressed in 3D polar co-ordinates:
$$mathbf{hat{n}} = mathbf{r}/|mathbf{r}| = (xmathbf{i} + ymathbf{j} + zmathbf{r} )/ |mathbf{r}| = sinphi costheta mathbf{i}+sinphi sintheta mathbf{j}+cosphi mathbf{k}$$
our logic above the included diagram should still hold: $dS = dA/mathbf{hat{n}}cdotmathbf{k}$ (here projecting onto the $xy$ i.e. $r,theta$ plane).
This yields $dS = dA/cosphi$ for our 3d polar case, which has got me stuck on the next step of evaluating
$$iint dS = iint secphi dA = iint sec(phi) rdrdtheta$$
It is intuitively true that as the angle $phi$ increases, our surface area element on the hemisphere becomes more and more vertical w.r.t. the $xy$ plane and so this $1/cosphi$ makes our corresponding area element larger, but how should I actually evaluate the integral?
$iint secphi dA = iint sec(phi) rdrdtheta$ is certainly ringing alarm bells, which makes me think that projecting on the $r,theta$ plane in the first place was a mistake, which is a problem not encountered when working in 3D cartesians.
Is there such a thing as projecting onto the “$theta,phi$ plane” and then integrating over those two variables? – if this were to be a valid method I believe it could solve my problem here. Please do ask for further clarification if it is needed.