If I take the upper hemisphere of a sphere, $x^2 + y^2 + z^2 = 1$, to be $\sqrt{1 – x^2 – y^2}$, then the normal is given by $\langle -f_x, – f_y, 1\rangle$ at any point.
This leads to an odd result: on the plane $z = 0$, while one might expect all normals of the sphere to not have any $z$ component (i.e, to only point radially outwards), the $\vec{k}$ component is still 1.
A similar parametrization in cylindrical coordinates is:
$\langle \sin(\phi)^2 \cos(\theta), \sin(\phi)^2 \sin(\theta), \sin(\phi) \cos(\phi)\rangle$.
At $\phi = 0$ and $\theta = 0$, which corresponds to the point $(0,0,1)$ in cartesian coordinates, the normal is $\langle 0,0,0 \rangle$ while one would expect $(0,0,1)$.
Best Answer
I think a large part of the difficulty here is the variation in the magnitudes of your "normals" in general.
When you give a "normal" in the form $\langle -f_x, - f_y, 1\rangle$, basically what you have is a radial projection of the hemisphere parameterized by $\langle \sin\phi \cos\theta, \sin\phi \sin\theta, \cos\phi\rangle$ (for $0 \leq \phi \leq \frac\pi2$) onto the plane $z = 1$. That is, you have a projection that takes a point $P$ on the surface of the hemisphere to a projected point $P'$ on the plane such that $P$ and $P'$ are collinear with the center of the sphere. (This is also called a gnomonic projection of the sphere.) The coordinates of your "normals" are the coordinates of the projected points.
As long as $0 \leq \phi < \frac\pi2$, each point of your sphere does in fact project onto that plane, although as $\phi$ approaches $\frac\pi2$ the magnitudes of your "normals" grow without bound. And obviously the boundary of your hemisphere, where $\phi = \frac\pi2$, does not project on to the plane $z = 1$ at all.
If you change your choice of $f$ as recommended in the answer by H.R., this problem goes away. All your normals then will have the same magnitude, and they will be defined at all points of the hemisphere.
When you give your "normals" in the form $\langle\sin^2\phi\cos\theta, \sin^2\phi\sin\theta, \sin\phi\cos\phi\rangle$, you again have non-uniform magnitudes, but this time the magnitudes go to zero as $\phi$ approaches zero. In effect, you are radially projecting the hemisphere onto a kind of degenerate torus given by $r = \sin\phi$.
But notice that the three components of your "normals" all have the common factor $\sin\phi$. You can normalize the magnitudes of all your "normals" (except for the case $\phi = 0$) by multiplying by the scalar $1/\sin\phi$. If you do this, you get the vectors $\langle\sin\phi\cos\theta, \sin\phi\sin\theta, \cos\phi\rangle$, that is, the coordinates of each vector (for $0 < \phi \leq \frac\pi2$) are simply the coordinates of the point on the sphere. If you define these vectors as the normals for all points such that $0 < \phi \leq \frac\pi2$, you may see that you can use the same formula to define the normal for $\phi = 0$ as well, and it works very nicely.
The formula ${\bf r}_\phi \times {\bf r}_\theta$ works nicely for the measurement of area (which is what it is used for in the example you quoted) precisely because its magnitude does go to zero as $\phi$ goes to zero and does so in the same way that the area element $r \,d\phi\,d\theta$ goes to zero. But as you observed, this property is not so desirable when you are trying to construct a set of normal vectors rather than trying to integrate some scalar function over area.