Yes, you can do this without any trigonometric functions, and even without any square roots. Take a look at Rodrigues' formula for a rotation matrix given an angle $\theta$ and an axis along a unit vector $k$:
$$R=I\cos\theta+\sin\theta[k]_\times+(1-\cos\theta)kk^{\text T}\;,$$
where $I$ is the identity matrix and
$$[k]_\times=\begin{pmatrix}0&-k_z&k_y\\k_z&0&-k_x\\-k_y&k_x&0\end{pmatrix}\;.$$
You have not an axis and an angle, but a vector $N$ that you want to rotate the $z$ axis into. The axis should be perpendicular to this vector and the $z$ axis, and the angle should be the angle between this vector and the $z$ axis. Since the $z$ component of your vector is the cosine of the angle it forms with the $z$ axis, you already have $\cos\theta=N_z$. Now consider
$$N_\perp=\begin{pmatrix}0\\0\\1\end{pmatrix}\times N=\begin{pmatrix}-N_y\\N_x\\0\end{pmatrix}\;.$$
This vector is perpendicular to both $e_z$ and $N$, so its direction is along the desired rotation axis. Its magnitude is $\sin\theta$. So
$$
\begin{eqnarray}
N_\perp&=&\sin\theta k\;,\\
[N_\perp]_\times&=&\sin\theta[k]_\times\;,\\
\end{eqnarray}
$$
where $k$ is a unit vector along the rotation axis. Now we have most of the ingredients of Rodrigues' formula. The only problem left is that we don't have $kk^{\text T}$ but $\sin^2\theta kk^{\text T}$. But $\sin^2\theta=1-\cos^2\theta$, and we have $\cos\theta$, so we can correct for that using only elementary arithmetic. Putting it all together, we have
$$
\begin{eqnarray}
R
&=&
I\cos\theta+[N_\perp]_\times+(1-\cos\theta)N_\perp N_\perp^{\text T}/(1-\cos^2\theta)
\\
&=&
I\cos\theta+[N_\perp]_\times+N_\perp N_\perp^{\text T}/(1+\cos\theta)
\\
&=&
IN_z+
\begin{pmatrix}0&0&N_x\\0&0&N_y\\-N_x&-N_y&0\end{pmatrix}
+\frac1{1+N_z}
\begin{pmatrix}N_y^2&-N_xN_y&0\\-N_xN_y&N_x^2&0\\0&0&0\end{pmatrix}
\\
&=&
\begin{pmatrix}N_z&0&N_x\\0&N_z&N_y\\-N_x&-N_y&N_z\end{pmatrix}
+\frac1{1+N_z}
\begin{pmatrix}N_y^2&-N_xN_y&0\\-N_xN_y&N_x^2&0\\0&0&0\end{pmatrix}
\;.
\end{eqnarray}
$$
You can check that this is an orthogonal matrix that rotates a unit vector along the $z$ axis into $N$. If you organize the operations efficiently, you only need two negations, three additions, one division and five multiplications to calculate the rotation matrix.
Note that the result becomes undefined at $N_z=-1$, corresponding to the fact that there's no preferred choice of axis in this case. You can either arbitrarily choose an axis
specifically for that case, e.g. the $x$ axis, or, since this might also cause numerical problems near $N_z=-1$, you can instead solve the problem for $-N_z$ whenever $N_z<0$, and then invert the resulting vectors. This will produce a mirror image of your distribution, but I'd assume that this is rotationally invariant around the $z$ axis anyway, and in that case the inversion wouldn't make a difference.
Note: one need to be careful when working in spherical coordinates, as there are two conventions used for the notation of the polar and azimuthal angles. See https://en.wikipedia.org/wiki/Spherical_coordinate_system and the 2 pictures near the top. In this answer I am using the physicist's convention, since it complies better with your own notation.
To clear the first point of confusion: All vectors start at the origin. They can be added together visually by placing them head to tail, but you should always think of them as having the tail at the origin and the head pointing in whatever direction it's supposed to.
The cross product in rectangular coordinates and the cross product in spherical coordinates are the same thing. The only difference is the way we represent them in formulas.
Notation also seems to be confusing you. The basis in spherical coordinates works differently from the basis in rectangular coordinates. See http://mathworld.wolfram.com/SphericalCoordinates.html and scroll down to the listing of unit vectors for a complete account. Caution: This article uses the mathematician's notation. You will have to swap the roles of $\theta$ and $\phi$ in your head, or just write things down to keep track. e.g. $\hat{\phi}$ there is what we'll think of as $\vec{e}_\theta$.
Here is a working definition of the cross product:
The cross product of two vectors $\vec{a}$ and $\vec{b}$ is another vector, mutually perpendicular to $\vec{a}$ and $\vec{b}$, with the direction determined by the right-hand rule. The magnitude of $\vec{a}\times\vec{b}$ is given by the area of the parallelogram spanned by $\vec{a}$ and $\vec{b}$.
Note that this definition makes no reference to the choice of coordinate system; it works in rectangular and spherical coordinates.
Therefore, if $\vec{a} = (a,\theta,\phi)$ and $\vec{b} = (b,0,0)$ in spherical coordinates, to find the cross product we need to find the direction the vector is pointing and the magnitude from this information. The key to this is to find the angle between $\vec{a}$ and $\vec{b}$. Well, $\vec{b}$ points straight along the $z$-axis, so it suffices to find the angle of $\vec{a}$ relative to the $z$-axis, and if you draw a picture you should see that this angle is $\theta$. (It doesn't matter what $\phi$ is here.) The area of the parallelogram spanned by $\vec{a}$ and $\vec{b}$ is then $ab\sin\theta$. As for the direction, to be perpendicular to the $z$-axis the vector must lie in the $xy$-plane, so we need to find the direction in the $xy$-plane perpendicular to $(a,\theta,\phi)$. Then it suffices to find the direction perpendicular to $(1,\pi/2,\phi)$, which is of course $(1,\pi/2,\phi\pm\pi/2)$. The right-hand rule says we take the minus sign, so
$$
\vec{a}\times\vec{b} = ab\sin\theta(1,\pi/2,\phi-\pi/2).
$$
We can take the plus sign by pointing in the opposite direction, so
$$
\vec{a}\times\vec{b} = -ab\sin\theta(1,\pi/2,\phi+\pi/2).
$$
Finally, let's express the last vector in rectangular coordinates. It is
$$
(1,\pi/2,\phi+\pi/2)_{\text{sph}} = (\sin(\pi/2)\cos(\phi+\pi/2),\sin(\phi+\pi/2)\sin(\pi/2), \cos(\pi/2))_{\text{rec}} = (-\sin\phi,\cos\phi,0)_{\text{rec}} = \vec{e}_{\phi}
$$
according to the definition of unit vectors in spherical coordinates. (See the mathworld link above if you haven't already.) Therefore we finally obtain
$$
\vec{a}\times\vec{b} = -ab\sin\theta\vec{e}_\phi
$$
as claimed.
Now to compute $\vec{a}\times\vec{a}\times\vec{b}$, we know that $\vec{a}$ and $\vec{a}\times\vec{b}$ are perpendicular, so they span a rectangle. Therefore the area they span is the product of their lengths, or $a^2b\sin\theta$. (And there is a minus sign from the previous computation.) For the direction we need to compute
$$
\frac{\vec{a}}{\|\vec{a}\|}\times\vec{e}_{\phi}.
$$
$\vec{e}_\phi$ is the vector lying in the $xy$-plane perpendicular to $(1,\pi/2,\phi)_{\text{rec}}$. To be perpendicular to this and $\vec{a}$, we must find the vector in the plane perpendicular to $\vec{e}_\phi$ which, if you draw a picture, you'll see is the one containing $\vec{a}$ and the $z$-axis. Therefore to find the perpendicular we just have to shift $\theta$ by $-\pi/2$ (the sign comes from the right-hand rule). Therefore
$$
\frac{\vec{a}}{\|\vec{a}\|}\times\vec{e}_{\phi} = (1,\theta-\pi/2,\phi)_{\text{sph}}
$$
and if we work this out in rectangular coordinates,
$$
(1,\theta-\pi/2,\phi)_{\text{sph}} = (\sin(\theta-\pi/2)\cos\phi,\sin(\theta-\pi/2)\sin\phi,\cos(\theta-\pi/2))_{\text{rec}} = (-\cos\theta\cos\phi,-\cos\theta\sin\phi,-\sin\theta)_{\text{rec}} = -\vec{e}_\theta.
$$
So putting everything together,
$$
\vec{a}\times\vec{a}\times\vec{b} = -a^2b\sin\theta(-\vec{e}_\theta) = a^2b\sin\theta\vec{e}_\theta.
$$
This is never zero except in the following cases (which are not mutually exclusive): if $\vec{a}=0$, if $\vec{b}=0$, or if $\vec{a}$ and $\vec{b}$ point in the same or opposite directions (so that $\sin\theta=0$). Which you could guess from the original working definition of the cross product: since $\vec{a}\times\vec{b}$ always points perpendicular to $\vec{a}$, as long as we don't have some trivial case the vectors $\vec{a}$ and $\vec{a}\times\vec{b}$ will span a parallelogram of nonzero area.
As you can see, this whole business with spherical coordinates is rather cumbersome for the cross product. This is because the cross product is designed to respect linear combinations, and therefore it works well in a rectangular system where the local geometry looks the same everywhere. However, in curvilinear coordinates (like spherical), vectors do not add coordinatewise like in rectangular coordinates, so the cross product cannot use the special linear structure of $\mathbb{R}^3$. The easiest way to do things is to learn to break up vectors into the local basis for spherical coordinates, and learn the formulas for the cross product of the local orthonormal basis in spherical coordinates. Either that, or as Griffiths suggests in Introduction to Electrodynamics, just convert to rectangular coordinates and try to do most things there.
Best Answer
Why not just reverse the ones with dot product <0 instead of discarding them? You will still have a uniform distribution in the hemisphere.