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?
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.
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.
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.
I was trying to understand this today and I appreciated the advice to refer to Khinchin's very clear analysis. I made a picture (above) that helped me see the key qualitative idea.
Consider a small volume of phase space trapped between two energy hypersurfaces, evolving with the Hamiltonian flow. Liouville's theorem tells us that the Lebesgue measure of the two volumes shown, $V_1 = \ell_1 w_1$ and $V_2 = \ell_2 w_2$, are equal. The $x_a$ coordinate in the picture represents all $6N-1$ dimensions of the energy surface $H^{-1}(E)$, while $x_b$ represents a single direction normal to the energy surface.
Since the dynamics of a time-independent Hamiltonian system with energy $E$ is constrained to the energy surface, we are interested in defining a measure on that $6N-1$-dimensional space that is invariant under the flow. But as can be seen in the schematic (blue vs. red strokes), the area on that surface is not invariant: $\ell_1 \ne \ell_2$.
Instead we can use the fact that, for small $\Delta E$, $w_i = \left.\frac{{\rm d} x_b}{{\rm d}E}\right|_i \,\Delta E$, and therefore: \begin{align} \ell_1 w_1 &= \ell_2 w_2\\ \ell_1 \left[ \left.\frac{{\rm d} x_b}{{\rm d}E}\right|_1 \,\Delta E\right] & = \ell_2 \left[\left.\frac{{\rm d} x_b}{{\rm d}E}\right|_2 \,\Delta E\right]\\ \frac{\ell_1}{\left. \frac{{\rm d} E}{{\rm d}x_b} \right|_1} &= \frac{\ell_2}{\left. \frac{{\rm d} E}{{\rm d}x_b} \right|_2}\\ \frac{\ell_1}{\left|\nabla H\right|_1} &= \frac{\ell_2}{\left|\nabla H\right|_2} \end{align} So, in this simplified schematic, phase space volume preservation due to Liouville implies preservation of $\int \frac{{\rm d}\sigma}{|\nabla H|}$.
