There are two ways of looking at this.
First, a rather rough "geometrical" look. Imagine an infinitesmal "rectangle" on the sphere having "top" at \phi and "bottom" at \phi+ d\phi, "left" at \theta and "right" at \theta+ d\theta. The left and right, "lines of longitude", are great circles, with radius the radius of the sphere, 1. The length of those two sides is 1(d\phi)= d\phi. The top and bottom, "lines of latitude", are not great circles. Their centers lie on the line through the poles, (0,0,1) and (0,0,-1). Given a point on one of those circles, drop a perpendicular to that vertical line. You have a right triangle with one leg being the radius of the circle, r, the hypotenuse the radius of the circle, 1, and angle at the center \phi. Then r/1= sin(\phi) so the radius of the circle is sin(\phi) and the length of the arc on the sphere is sin(\phi)d\theta. The area of the "rectangle" is (width times height) sin(\phi)d\theta d\phi.
A more rigorous calculation involves the "fundamental vector product" which is worth knowing on its own. Any smooth surface, since a surface is two dimensional, can be written in terms of two parameters: x= f(u,v), y= g(u,v), z= h(u,v) which can also be written as the vector equation: \vec{r}= f(u,v)\vec{i}+ g(u,v)\vec{j}+ h(u,v)\vec{k}. The derivatives of that with respect to the two variables, \vec{r_u}= f_u(u,v)\vec{i}+ g_u(u,v)\vec{j}+ h_u(u,v)\vec{k} and \vec{r_v}= f_v(u,v)\vec{i}+ g_v(u,v)\vec{j}+ h_v(u,v)\vec{k} are vectors in the tangent plane to the surface at each point. The cross product then, \vec{r_u}\times\vec{r_v}, the "fundamental vector product", is perpendicular to the surface and, because of the derivatives, its length measures the area of an infinitesmal region: dS= \left|\vec{r_u}\times\vec{r_v}\right|dudv.
The spherical coordinates are connected to the cartesian coordinates by x= \rho cos(\theta)sin(\phi), y= \rho sin(\theta)sin(\phi), z= \rho cos(\phi) and we can take those as parametric coordinates for the unit sphere by taking \rho= 1: x= cos(\theta)sin(\phi), y= sin(\theta)sin(\phi), z= cos(\phi) or \vec{r}= cos(\theta)sin(\phi)\vec{i}+ sin(\theta)sin(\phi)\vec{j}+ cos(\phi)\vec{k}. The derivatives are \vec{r_\theta}= -sin(\theta)sin(\phi)\vec{i}+ cos(\theta)sin(\phi)\vec{j} and \vec{r_\phi}= cos(\theta)cos(\phi)\vec{i}+ sin(\theta)cos(\phi)\vec{j}- sin(\phi)\vec{k}. The cross product of those is cos(\theta)sin^2(\phi)\vec{i}+ sin(\theta)sin^2(\phi)\vec{j}+ sin(\phi)cos(\phi)\vec{k} which has length \sqrt{cos^2(\theta)sin^4(\phi)+ sin^2(\theta)sin^4(\phi)+ sin^2(\phi)cos^2(\phi)}= \sqrt{sin^4(\phi)+ sin^2(\phi)cos^2(\phi)}= \sqrt{sin^2(\phi)(sin^2(\phi)+ cos^2(\phi))}= sin(\phi). That gives dS= sin(\phi)d\phi d\theta.
By the way, my \theta and \phi are reversed from yours. Mathematics notation (mine) uses \theta as "longitude" and \phi as "co-latitude". Engineering notation reverses those.
