What's the intuition behind the Co-Area formula?
This is most easily understood when $N = 1$. Suppose $A \subset \mathbb{R}^m$ is a bounded domain and $f: A \rightarrow \mathbb{R}$ is a smooth function whose gradient is everywhere nonzero.
Its level sets are non-intersecting hypersurfaces. Suppose first that you want to compute the volume of the $A$ in terms of the surface areas of the hypersurfaces. Suppose $a$ is the inf of $f$ on $A$ and $b$ is the sup. You can divide the interval $[a,b]$ into equal sized subintervals $I_1 = [t_0,t_1], \dots, I_N = [t_{N-1},t_N]$ of size $\delta = (b-a)/N$. The volume of $A$ is the sum of the volumes of $f^{-1}(I_k)$. On the other hand, each $f^{-1}(I_k)$ is a shell with varying thickness. At each point, the thickness is roughly $\Delta t/|\nabla f|$. So the volume of the shell is roughly $$ V(f^{-1}(I_k)) \simeq \Delta t\int_{f^{-1}(t_k)}\frac{dA}{|\nabla f|}, $$ where $dA$ is the $(M-1)$-dimensional Hausdorff measure on $f^{-1}(t_k)$. Adding up the volumes of the shells and taking the limit $N \rightarrow \infty$, we see that the volume of $A$ is $$ V(A) = \int_{-\infty}^{\infty} \int_{f^{-1}(t)} \frac{dA_t}{|\nabla f|}\,dt, $$ where $dA_t$ is the $(M-1)$-dimensional Hausdorff measure of $f^{-1}(t)$.
If, on the other hand, you integrate $|\nabla f|$ on $A$ using the same approach, you get $$ \int_{f^{-1}(I_k)} |\nabla f|\,dx \simeq \Delta t\int_{f^{-1}(t_k)}|\nabla f|\frac{dA}{|\nabla f|}, $$ and therefore $$ \int_A |\nabla f|\,dx = \int_{-\infty}^{\infty} \int_{f^{-1}(t)}\,dA_t\,dt $$ The $N>1$ case is similar, except you chop $\mathbb{R}^N$ into small rectangular pieces and use the fact that the cross section of the inverse image of each piece is roughly a parallelogram whose volume is $(J_Nf)^{-1}$. and therefore $$ \int_A J_Nf\,dx = \int_{\mathbb{R}^N} \int_{f^{-1}(y)}\,dA_y\,dy, $$ where $dA_y$ is the $(M-N)$-dimensional Hausdorff measure on $f^{-1}(y)$.
The general formula is $$ \int_A \phi(x)\,dx = \int_{\mathbb{R}^N} \int_{f^{-1}(y)} \phi(x)\frac{dA_y}{J_Nf}\,dy $$