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}}\cdot\mathbf{k}\, dS \rightarrow dS = dA/\mathbf{\hat{n}}\cdot\mathbf{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 2\pi, 0 \leq r \leq a $) to reach an answer of $2\pi 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}\| = (x\mathbf{i} + y\mathbf{j} + z\mathbf{r} )/ \|\mathbf{r}\| = \sin\phi \cos\theta \mathbf{i}+\sin\phi \sin\theta \mathbf{j}+\cos\phi \mathbf{k}$$
our logic above the included diagram should still hold: $dS = dA/\mathbf{\hat{n}}\cdot\mathbf{k}$ (here projecting onto the $xy$ i.e. $r,\theta$ plane). This yields $dS = dA/\cos\phi$ for our 3d polar case, which has got me stuck on the next step of evaluating
$$\iint dS = \iint \sec\phi dA = \iint \sec(\phi) rdrd\theta$$
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/\cos\phi$ makes our corresponding area element larger, but how should I actually evaluate the integral?
$\iint \sec\phi dA = \iint \sec(\phi) rdrd\theta$ 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.
First of all, they are usually called spherical coordinates, not polar, when you are in $3D$. The radial variable is usually called $\rho$ instead of $r$. The element of area on the surface of a sphere is the product of two infinitesimal lengths. $\rho d\phi$ is the "changing latitude" length, and we usually say $r = \rho \sin \phi$ for the horizontal radius at that "latitude", so that $\rho \sin \phi d\theta$ is the "east-west" length, giving a surface area element $dS = \rho^2\sin\phi d\phi d\theta$. On this surface, $\rho = a$. So yes, you can absolutely integrate that way.