TL;DR: The logic behind the equation $|\psi\rangle=\cos\frac{\theta}{2}|0\rangle+e^{i\phi}\sin\frac{\theta}{2}|1\rangle$ can be iterated to obtain a real parameterization for states of arbitrary finite dimension.
Iterative procedure
We parameterize an expansion of $|\psi\rangle$ in basis $|0\rangle,\dots,|N\rangle$ in $N$ iterations. The $k$th iteration for $k=1,\dots,N$ consists of two steps. First, we introduce $\theta_k\in[0,\pi]$ to apportion the unit norm between the magnitudes of the amplitudes of $|k-1\rangle$ and the projection $|\psi_{k\dots N}\rangle$ of $|\psi\rangle$ onto the subspace spanned by $|k\rangle,\dots,|N\rangle$. Next, we stick a relative phase $\alpha_k\in[0,2\pi)$ on $|k\rangle$. Thus, the $k$th iteration rewrites $|\psi_{k-1\dots N}\rangle$ as
$$
|\psi_{k-1\dots N}\rangle=\cos\frac{\theta_k}{2}|k-1\rangle + e^{i\alpha_k}\sin\frac{\theta_k}{2}|\psi_{k\dots N}\rangle.\tag1
$$
The final result can be simplified by setting $\phi_k:=\alpha_1+\dots+\alpha_k\mod 2\pi$.
Each step introduces a single real parameter, for a total of $2N$ real parameters, as expected for an $N+1$ dimensional system.
Two-qubit case
For a two-qubit state the above procedure yields
$$
\begin{align}
|\psi\rangle&=\cos\frac{\theta_1}{2}|00\rangle\\
&+ e^{i\phi_1}\sin\frac{\theta_1}{2}\cos\frac{\theta_2}{2}|01\rangle\\
&+ e^{i\phi_2}\sin\frac{\theta_1}{2}\sin\frac{\theta_2}{2}\cos\frac{\theta_3}{2}|10\rangle\\
&+ e^{i\phi_3}\sin\frac{\theta_1}{2}\sin\frac{\theta_2}{2}\sin\frac{\theta_3}{2}|11\rangle
\end{align}\tag2
$$
where $\theta_1,\theta_2,\theta_3\in[0,\pi]$ and $\phi_1,\phi_2,\phi_3\in[0,2\pi)$.
This isn't the only way to parameterize a two-qubit state. An alternative approach uses Schmidt decomposition and single-qubit parameterization.