Let $S$ be a surface given by the parametrization $\vec{r}(u,v)= x(u,v)\vec{i} + y(u,v)\vec{j} + z(u,v)\vec{k}$ for $(u,v)\in D$. If $f(x,y,z)$ is a function defined on $S$, then the surface integral can be found using the formula $$ \iint_S f(x,y,z) dS = \iint_D f(\vec{r}(u,v))\lvert \vec{r}_u\times\vec{r}_u\rvert dA $$
I have seen how to derive this by drawing some pictures, but I am wondering if there is a proof that makes use of the change of variable formula. I have tried a few things, but I can't find anything that gives a Jacobian of $\lvert \vec{r}_u\times\vec{r}_u\rvert$.
Is this possible?
The formula you have written down is merely the definition of the surface integral(and as such doesn't need any proof). Ok, to be slightly more precise, if you define the integral as the sum of surface "area elements" times the value of the function at any point inside that surface, then you see that this gives you the above formula after noting that the vector product is the same as the area of the surface element.
The change of variable formula just shows that you may parametrize your surface S in any way you like and still get the same answer. The Jacobian comes up when we have an integral of the type $\int_{a}^{b} \int_{c}^{d}f(x,y) \,dx\,dy$ and some other params $u,v$ and $u=u(x,y)$ and $v=v(x,y)$ and then $\int_{a}^{b} \int_{c}^{d}f(x,y) \,dx\,dy = \int_{r}^{s} \int_{m}^{n}J(u,v)f(u,v) \,dx\,dy$.
In our case, assuming you have some other set of parameters $(p,q)$ such that $\vec{r}(p,q) = x(p,q)\vec{i}+y(p,q)\vec{j}+z(p,q)\vec{k}$ , then $|\vec{r}_u \wedge \vec{r}_v|dudv= J(u,v)|\vec{r}_p \wedge \vec{r}_q|dpdq$(write down this computation in full if you don't believe me!) and by the change of variable formula, the two integrals(one in which we parametrized using $p,q$ and the other with $u,v$) will be equal.