Suppose I transform an integral $$I=\int f(x,y) \, dx \, dy$$ using polar coordinates, setting $x=r\cos\theta$ and $y=r\sin\theta$. We get $$ \begin{split} dx &= \cos\theta \, dr - r\sin\theta \, d\theta\\ dy &= dr\,\sin\theta + r \cos \theta \, d\theta. \end{split} $$
Now the volume element $dx \, dy$, should become $$ (\cos\theta \, dr - r\sin\theta \, d\theta)(dr\,\sin\theta + r \cos \theta \, d\theta), $$ whereas my book says it becomes $r \, dr \, d\theta$.
I understand the geometric argument but I am looking for a way to connect the result the result of my book with what should I have got had I used differentials.
Multivariable transformations don't quite work out so easily in multi-d to a simple 1d change. You have to generalize the 1D change of variables carefully (see here for example).
You really have the 2D transformation of coordinates $$ \begin{pmatrix} dx \\ dy \end{pmatrix} = \begin{pmatrix} \cos \theta & -r \sin \theta \\ \sin \theta & r \cos \theta \end{pmatrix} \begin{pmatrix} dr \\ d\theta \end{pmatrix} $$ so the correct change of variable factor is given by the Jacobian, which is the determinant of the transformation matrix. You get $$ \det \begin{pmatrix} \cos \theta & -r \sin \theta \\ \sin \theta & r \cos \theta \end{pmatrix} = r \cos^2 \theta + r \sin^2 \theta = r. $$