Indeed, the distribution function of the eigenphases of a random matrix in $\operatorname{SO}(n)$ has a peak at 0 and at $\pm\pi$. It only becomes uniform for large $n$. The joint distribution function for $n=2m$ even is given by (Girko, 1985) $$p(\theta_1, \cdots, \theta_m) = C \prod_{1 \leq k < j \leq m} (\cos\theta_k - \cos\theta_j)^2~.$$ As a simple example, for $n=4$ this gives upon integration over $\theta_2$ at fixed $\theta_1$ the density profile $$p(\theta)=\frac{\cos 2 \theta+2}{4 \pi }.$$ If you sample from $O(n)$, you can either discard the matrices with determinant $-1$, or multiply the first column by $-1$.
By way of illustration, this is the density profile for $n=100$ (computed by averaging over 500 random matrices), when the peaks at $0$ and $\pm\pi$ have become very small.
The reason why the Haar measure for SO$(n)$ does not give a rotationally invariant density of eigenvalues on the unit circle is that the eigenvalues must come in complex conjugate pairs $e^{\pm i\theta}$. This constraint is not rotationally invariant.