I'll spare you the cotangent bundles and differential geometrese and just summarize the takeaway that, in fact, $J=1/r^2 \sin \theta$ so that the entire phase-space Jacobian is 1,
$$
dx dy dz ~ dp_x dp_y dp_z = dr d\theta d\phi ~dp_r dp_\theta dp_\phi .
$$
A direct (blood, sweat and tears) derivation is available in Peter Joot's Blog.
The reason is that Cartesian to spherical is a point canonical transformation, so it preserves phase-space volumes (Liouville's theorem―which also holds for motion, since that is also a canonical transformation generated by the Hamiltonian).
To rationalize this, consider a free particle of mass m =1. The Hamiltonian is then $\vec p ^2/2$, generating
$$
\frac{d\vec r}{dt} = \{\vec r , \vec p ^2 \}/2 = \vec p.
$$
It is simple in Cartesian coordinates, but in spherical coordinates,
given the line element
$$
d\vec r= \hat r ~ dr +\hat \theta ~ r d\theta + \hat \phi ~ r \sin\theta d\phi,
$$
you have
$$
\dot{\vec r}= \hat r ~ \dot{r} +\hat \theta ~ r \dot{\theta} + \hat \phi ~ r \sin\theta ~\dot{\phi} \\
=\vec p= \hat r ~ p_r +\hat \theta ~ \frac{1}{r} p_\theta + \hat \phi ~\frac{1}{ r \sin\theta} p_\phi ~~.
$$
These are the canonical conjugate momenta gotten from the canonical procedure and, e.g., $p_\phi$ is not the projection of $\vec p$ in the direction $\phi$!
You've seen this covariant bit before in the gradient expressed in spherical coordinates,
$$
\nabla = \hat r ~ \partial_r +\hat \theta ~ \frac{1}{r} \partial_\theta + \hat \phi ~\frac{1}{ r \sin\theta} \partial_\phi ,
$$
not coincidentally, as it is proportional to the quantization of the momentum when you transcend classical mechanics.
Twice the Hamiltonian is, in this language,
$$
2H= \vec p^2 = p_r^2 + \frac{p_\theta^2}{r^2} + \frac{p_\phi^2}{r^2\sin^2\theta },
$$
(and the Liouville one-form would be $\vec p \cdot d\vec r = p_r dr + p_\theta d\theta +p_\phi d\phi $. In components, $p_r=\dot{r}, \quad p_\theta /r = r \dot{\theta}, \quad p_\phi/ r \sin\theta= r\sin \theta ~ \dot{\phi} $. )
The volume element in phase space, then, by above, is
$$
d^3 \vec r ~ d^3 \vec p= r^2 \sin \theta ~ dr d\theta d\phi ~ \frac{1}{r^2 \sin \theta} dp_r dp_\theta dp_\phi,
$$
collapsing to the top line.
