Use 'polar' coordinates, as in $\phi(\lambda, \theta) = (\lambda a \cos \theta, \lambda b \sin \theta)$, with $(\lambda, \theta) \in S = (0,1] \times [0,2 \pi]$. It is straightforward to compute the Jacobian determinant as
$$ J_{\phi}(\lambda, \theta) = |\det D\phi(\lambda, \theta)| = \lambda a b.$$
Let $E = \{ (x,y) \,|\, 0 <(\frac{x}{a})^2 + (\frac{y}{b})^2 \leq 1 \}$. (Eliminating $(0,0)$ makes no difference to the integral, and is a technicality for the change of variables below.) We have $E = \phi (S)$, and
$$\begin{align}
I &= \rho \int_{\phi ( S)} (x^2+y^2) \, dx dy \\
&= \rho \int_{S} \lambda^2 (a^2 \cos^2 \theta+ b^2 \sin^2 \theta) \lambda a b \, d \lambda d \theta \\
&= \rho a b \int_{0}^1 \lambda^3 \, d \lambda \int_0^{2 \pi} a^2 \cos^2 \theta+ b^2 \sin^2 \theta \, d\theta \\
&= \rho \pi a b \frac{a^2+b^2}{4}.
\end{align}$$
Best Answer
You are almost right. Note that the density $\rho$ is just proportional to $xy$, so we can write $\rho=Cxy$. To find the value of $C$, you know that the total mass is $M$: $$M=\int_0^a dx\int_0^\sqrt{b^2(1-\frac{x^2}{a^2})} Cxy dy$$ So $$C=\frac{M}{\int_0^a dx\int_0^\sqrt{b^2(1-\frac{x^2}{a^2})} xy dy}$$ Then your moment of inertia $I$ is: $$I=M\frac{\int_0^a dx\int_0^\sqrt{b^2(1-\frac{x^2}{a^2})} xy(x^2+y^2) dy}{\int_0^a dx\int_0^\sqrt{b^2(1-\frac{x^2}{a^2})} xy dy}$$