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 (<A HREF="https://doi.org/10.1007%2Fbf01061167">Girko, 1985</A>)
$$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.

<IMG SRC="https://i.sstatic.net/yOs9r.png" WIDTH="300"/>

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.