The invariant measure on an energy surface of a Hamiltonian system Consider a Hamiltonian system with a time-independent Hamiltonian $H (p, q )$. By the Liouville theorem, the measure $d^np d^nq $ is conserved. 
However, one should also notice that the energy is conserved and the system does not evolve in a space, but on a hypersurface, i.e., the energy surface $E = H(p, q)$. 
So, what is the invariant measure on the energy surface, if there exists such a measure at all?   
 A: You are right that while Liouville's theorem says something about invariance of the measure $\Pi dq dp$, it does not say anything directly about an invariant measure on the hypersurface given by $E=H(p,q)$.  
In his excellent book Mathematical Foundations of Statistical Mechanics, Khinchin (Section 7) shows that there indeed exists an invariant measure on that hypersurface.  That invariant measure $\mu$ is given by:
\begin{equation}
\mu(M) = \int_M \frac{ds}{|\nabla H|},
\end{equation}
for a set $M$, where $ds$ is a volume element in that hypersurface.  As pointed out here roughly speaking, $1/|\nabla H|$ describes the thickness of the energy hypersurface at that point.
Added: And here is another thread on this topic.
A: As Menachem said, 
\begin{equation}
\mu(M) = \int_M \frac{ds}{|\nabla H|},
\end{equation}
is a good answer. 
Now I want to give a more abstract answer $\mu_E$ with surface element $d\sigma$ that will be described below.
Consider N particles with phase space measure $d\tau=dq_1...dq_N dp_1...dp_N$
Let $d\sigma$ be a measure on energy surface such that $d\sigma dE=d\tau$. Call the measure $\mu_E$, so $\mu_E(V)=\int_{V\cap \Sigma_E} d\sigma$. Here, $\Sigma_E$ is the energy surface and V is any phase space volumn. We can then prove such surface measure is hamiltonian-flow-invariant. 
Let $\phi_t()$ be the hamiltonian flow ( ie, a time evolution of initial state $\tau_0$ is $\phi_t(\tau_0)$ ). By Liouvill's Theorem, $d\tau$ is hamiltonian-flow-invariant. 
$$ \int_V d\tau =\int_{\phi_t(V)} d\tau $$
Now cut $V$ be between $E$ and $E+\Delta E$ and only keep the shell.
$$ \int_{E}^{E+\Delta E} \int_{V \cap \Sigma_{E'}} dE'd\sigma = \int_{E}^{E+\Delta E} \int_{\phi_t(V) \cap \Sigma_{E'}} dE'd\sigma$$ 
This holds true for any $E$ and $\Delta E$ that we choose (because of conservation of energy). Now we substitute the definition of $\mu_E(V)$ and $\mu_E(\phi_t(V))$ 
$$ \int_{E}^{E+\Delta E} dE' \mu_{E'}(V) = \int_{E}^{E+\Delta E} dE' \mu_{E'}(\phi_t(V))$$
Since this holds true for any $E$ and $\Delta E$, we divide both sides by $1/\Delta E$, fix E, then take $\Delta E$ to approach zero.
We get
$$ \mu_{E}(V)=\mu_E(\phi_t(V)) $$
The above derivation shows if we can construct a $d\sigma$ such that $d\sigma dE=d\tau$, then the construction is flow-invariant.
Now we show the measure Menachem wrote, $(ds/ |\nabla H| )$, does satisfy $(ds/ |\nabla H|) \ dE= d\tau $. Since $ds\ dh=d\tau$, where $dh$ is the height of the volume element in phase space, we want to show $dh=dE/ |\nabla H|$ or $dh |\nabla H|=dE$, and this is almost the definition of derivative. (If we consider $\vec{dh}$ as a vector on the phase space normal to the energy surface, then $dE=\vec{\nabla} H \cdot \vec{dh}=|\nabla H| dh$ ) 
In the forward direction, $dh |\nabla H|=dE$ -> $dE/|\nabla H|=dh$ -> $dE/|\nabla H|\ ds=dh\ ds=d\tau$ -> $(ds/|\nabla H|)\ dE=d\tau$. By the previous derivation $ds/|\nabla H|$ gives a flow invariant measure.
A: The Liouville measure is basically the volume of a tiny full-dimensional region in phase space.  If you know the position and momentum of every particle exactly, the the state of the system is described by a point in phase space, and Liouville's theorem becomes trivial, because obviously a point is always zero-dimensional no matter how you move it.  It's true that the point's trajectory will be confined to a constant-H hypersurface (if the Hamiltonian doesn't depend explicitly on time), but there's no interesting "measure" to be conserved.
Liouville's theorem becomes useful when you have a statistical ensemble of systems, or more prosaically, when the positions and momenta have some uncertainty.  For example, consider the ensemble of one-particle systems where the position and momenta are initially uniformly distributed over the intervals $(x, p) \in [x_0 - \epsilon, x_0 + \epsilon] \times [p_0 - \delta, p_0 + \delta]$, for some small $\epsilon$ and $\delta$.  Then the ensemble is described by an $\epsilon \times \delta$ rectangle in phase space, and the area of this region will distort over time but its area $\epsilon \delta$ will remain constant.  Such an ensemble will generically extend over a finite range of values $H$ of the Hamiltonian.  So in the ensemble context where Liouville's theorem is useful, "the system" is not in fact confined to a single hypersurface.
In a given coordinate system on phase space, you can define a codimension-1 hypersurface $A(E)$ by the intersection of the volume of the ensemble system and the hypersurface with energy $E$, and Liouville's theorem guarantees the conservation of the volume $\int_{-\infty}^\infty A(E)\, dE$.  But I'm not sure if the function $A(E)$ itself is conserved, because as you suggest, it's not obvious that this hypersurface area measure can be defined in a way that is left invariant under canonical transformations.
