To integrate this region in polar coordinates,
it is advisable to break up the integral into two parts,
as shown in the figures below:
The two parts of the integral are divided by the diagonal line through
the upper right corner of the rectangle.
Since the sides of the rectangle are $a$ and $b$,
this diagonal line is at the angle $\arctan \frac ba.$
For $0 \leq \theta \leq \arctan \frac ba,$
you would integrate over $0 \leq r \leq a \sec\theta,$
and for $\arctan \frac ba \leq \theta \leq \frac\pi2,$
you would integrate over $0 \leq r \leq b \csc\theta.$
If you actually try this, I think you'll find that it is no easier than
doing the integration in rectangular coordinates.
It may even be worse.
An alternative approach, rather than combining $x^2+y^2$ into $r^2$,
is to integrate the terms separately:
$$\begin{eqnarray}
m &=& \int_0^a\int_0^b (1+x^2+y^2)\,dy\,dx \\
&=& \int_0^a\int_0^b dy\,dx
+\int_0^a\int_0^b x^2 \,dy\,dx
+\int_0^a\int_0^b y^2 \,dy\,dx
\end{eqnarray}$$
$$\begin{eqnarray}
m\bar{x} &=& \int_0^a\int_0^b x(1+x^2+y^2)\,dy\,dx \\
&=& \int_0^a\int_0^b x \,dy\,dx
+\int_0^a\int_0^b x^3 \,dy\,dx
+\int_0^a\int_0^b xy^2 \,dy\,dx
\end{eqnarray}$$
$$\begin{eqnarray}
m\bar{y} &=& \int_0^a\int_0^b y(1+x^2+y^2)\,dy\,dx \\
&=& \int_0^a\int_0^b y \,dy\,dx
+\int_0^a\int_0^b x^2y \,dy\,dx
+\int_0^a\int_0^b y^3 \,dy\,dx
\end{eqnarray}$$
Now you have nine integrals to solve, but they're all quite simple.
Center of mass calculation is essentially the same as other averages such as expectation value in probability theory. Essentially it is the weighted average of the space coordinate.
The definition is best wrote down in Cartesian coordinates:
$$\mathbf{r}_\text{center} = \int_V \mathbf{r} f(\mathbf{r}) \mathrm{d} V$$
with $\mathbf{r} = (x,y,z)$ and $\mathrm{d} V = \mathrm{d} x \mathrm{d} y \mathrm{d} z$. And $$f(\mathbf{r}) = \frac{1}{\int_V\rho(\mathbf{r}) \mathrm{d} V} \rho(\mathbf{r})$$
the volume normalized distribution function.
This hint should allow you to solve the stated problem. It is now just a matter of correct calculation of the integrals with spherical coordinates. You will get the center of mass in Cartesian coordinates, though. But I dont if there exists any simple integral expression to obtain the center of mass directly in spherical coordinates.
Best Answer
You don't need to use Polar Coordinates.
With the same set-up as you have described, consider semi-circular elements of radius $x$ and thickness $\delta x$
Then the area element is $\delta A=\pi x\delta x$
Now if the density per unit area is $\rho =kx$, then the mass element is $$\delta m=\pi kx^2\delta x$$
Integrating gives the mass of the lamina, as you have found, is $$m=\frac 13 \pi kr^3$$
Now using the formula for the centroid of an arc, the distance of the centroid of the semicircular element from the origin is $$\frac {x\sin \left(\frac {\pi}{2}\right)}{\frac {\pi}{2}}=\frac {2x}{\pi}$$
Now applying Varignon's Principle, $$\frac 13 \pi kr^3\bar x=\int_0^r\frac{2x}{\pi}\pi kx^2dx=\frac{kr^4}{2}$$
Hence the result $$\bar x=\frac{3r}{2\pi}$$