You want to think of the Haar measure $d\mu(U)$ as a way of measuring uniformity in the group $U(N)$ of unitary $N\times N$ matrices.
To form your intuition, consider $N=1$. You then have $U=e^{i\phi}$, with $0<\phi\leq 2\pi$ and $d\mu(U)=d\phi$ measures the perimeter of the unit circle. This is a uniform measure, because $d(\phi+\phi_0)=d\phi$ for any fixed phase shift $\phi_0$. You could write the requirement of uniformity in the form $d\mu(UU_0)=d\mu(U)$, with $U_0=e^{i\phi_0}$ the unitary matrix corresponding to the phase shift $\phi_0$.
Once your intuition is formed for $N=1$, you simply generalize to $N>1$ using the same definition of uniformity, $d\mu(UU_0)=d\mu(U)$ for any fixed $U_0\in U(N)$. For orthogonal (or symplectic) matrices you use the same definition of uniformity, with $U_0$ now restricted to the orthogonal or symplectic subgroup of $U(N).$
To explicitly write down the Haar measure $d\mu(U)$ in terms of the matrix elements of $U$ is only easily done for a few small values of $N$. (In particular, there is no relationship to random directions of rows or columns, as Yemon Choi pointed out.) You typically do not need such explicit expressions, since integrals with the Haar measure can be evaluated by using only the definition of uniformity.
In response to the follow-up question: If you wish to evaluate Haar-measure integrals of polynomials of matrix elements of $U$, you can use the socalled Weingarten functions.
http://en.wikipedia.org/wiki/Weingarten_function
Here is a Mathematica program to generate these,
http://arxiv.org/abs/1109.4244
If you need an explicit expression for the Haar measure, the steps to take are the following:
1) parameterize your matrix $U$ in terms of a set of real parameters $\{x_i\}$.
2) calculate the metric tensor $m_{ij}$, defined by $\sum_{ij}|dU_{ij}|^2 = \sum_{ij}m_{ij}dx_i dx_j$
3) obtain the Haar measure by equating $d\mu(U) = ($Det $m)^{1/2}\prod_i dx_i$
This is the general recipe. In practice, for many parameterizations the answer is in the literature. In particular, for the Haar measure in Euler angle parameterizations see:
http://arxiv.org/abs/math-ph/0205016
http://www.cft.edu.pl/~karol/pdf/ZK94.pdf
Your general determinant has to do with the Wronskian isomorphism for $SL_2$ representations
namely a map
$$
\wedge^m({\rm Sym}^n(\mathbb{C}^2))\rightarrow {\rm Sym}^m({\rm Sym}^{n-m+1}(\mathbb{C}^2)))\ .
$$
If you follow it by the natural map
$$
{\rm Sym}^m({\rm Sym}^{n-m+1}(\mathbb{C}^2))\rightarrow {\rm Sym}^{m(n-m+1)}(\mathbb{C}^2)
$$
which corresponds to restricting to the diagonal, then you get the usual Wronskian or rather a homogenized version of it where polynomials in one variable are seen as binary forms.
You can find more details in this paper and that one.
Great question by the way and thank you for pointing out the article by Ping Sun.
Best Answer
The Poisson kernel for the orthogonal group was calculated by Benjamin Béri in Generalization of the Poisson kernel to the superconducting random-matrix ensembles. This is in the context where $X$ is the scattering matrix of a superconductor in the Cartan symmetry class D (broken time-reversal and spin-rotation symmetries).
The result for $X\in{\rm SO}(N)$ (real unitary matrices with determinant $+1$) is $$P_{\rm SO}(X)\propto\frac{1}{|\operatorname{det}(X-Z)|^{N-1}},$$ to be contrasted with $$P_{\rm U}(X)\propto\frac{1}{|\operatorname{det}(X-Z)|^{2N}}$$ for $X\in{\rm U}(N)$.
Since the exponent $N-1$ differs from the expectation in the OP, let me check the reproducing property for ${\rm SO}(2)$ and $Z=zI_2$ (with $-1<z<1$). In that case the Poisson kernel is $$P_{\rm SO}(X)d\mu(X)=f(\theta)\,d\theta,\;\;f(\theta)=\frac{1-z^2}{2 \pi \left(z^2-2 z \cos \theta+1\right)}, $$ $$\text{in the parameterization}\;\;X(\theta)=\begin{pmatrix}\cos\theta&\sin\theta\\ -\sin\theta&\cos\theta\end{pmatrix}\;\;0<\theta<2\pi,$$ and one readily checks that $$\langle X^p\rangle\equiv\int_0^{2\pi} X^p(\theta)f(\theta)\,d\theta=z^p I_2=Z^p.$$ The reproducing property for the orthogonal Poisson kernel fails if $Z$ is not proportional to the unit matrix, for example, if $Z={{z\;0}\choose{0\,-z}}$ one has $\operatorname{det}(X-Z)=1-z^2$, independent of $X$ and $\langle X\rangle=0\neq Z$.
To generalize the Poisson kernel to all three classical groups ${\rm U}(N)$, ${\rm SO}(N)$, ${\rm Sp}(N)$, one defines an $N\times N$ subunitary matrix $Z$ as the submatrix of the $2N\times 2N$ unitary matrix $$\Omega=\begin{pmatrix}Z&T'\\ T&Z'\end{pmatrix}.$$ The $N\times N$ unitary matrix $X$ is then constructed by $$X=Z+T'X_0(1-Z'X_0)^{-1}T,$$ with $X_0$ an $N\times N$ unitary matrix distributed according to the Haar measure on the unitary, orthogonal, or symplectic group. The Poisson kernel is then defined as the corresponding probability distribution function $P(X)$ for a given $Z$. The result is $P(X)\propto|\operatorname{det}(X-Z)|^{-\beta N+\gamma}$, with $\beta=2,\gamma=0$ for ${\rm U}(N)$, $\beta=1,\gamma=1$ for ${\rm SO}(N)$, and $\beta=1,\gamma=-1$ for ${\rm Sp}(N)$. $$\mbox{}$$ The reproducing property, $\int P(X)f(X)d\mu(X)=f(Z)$, is ensured if the average of the matrix product $X_0(Z'X_0)^p$ vanishes for all $p=0,1,2\ldots$, which happens for any integer $N$ and any $Z$ in the unitary group, and for even $N$ and $Z\propto I_N$ in the orthogonal or symplectic groups.