The reason why cylindrical coordiates are natural for this problem is because, as you have noted with your graphed picture of the domain, you are dealing with part of a cylinder. The key thing to notice here, however, is which axis is the axis of rotation. Now I couldn't clearly tell what the region you're saying we are integrating over is, as I think you wrote it in a rather unclear way, but I assume the region we are integrating over is the set
$$\Omega = \{(x,y,z)\in\mathbb{R}^3 : x^2+z^2\leq1 \land 0\leq y\leq 2 \land z \geq 0\}.$$
If this is the case, then you can see from the part $x^2+z^2\leq 1$ that the axis of rotation for the cylinder is the $y$-axis. Therefore it is natural to pick the cylindrical coordinates
$$\begin{cases}
x=r\cos \theta, \\
y=y,\\
z=r\sin \theta,
\end{cases}$$
as this corresponds to a cylinder with the $y$-axis as its axis of rotation. We now investigate what $\Omega$ corresponds to over this coordinate transform. So we consider the definition of $\Omega$ piece by piece. First we have that
$$x^2+z^2\leq 1 \iff r^2 \cos^2\theta +r^2\sin^2\theta=r^2(\cos^2\theta+\sin^2\theta)=r^2\leq 1,$$
so $0\leq r \leq 1$. Secondly we have that $0\leq y\leq 2$. Lastly we have that (with $r\neq 0$ and as $0\leq \theta\leq 2\pi$)
$$z\geq 0 \iff r\sin \theta\geq 0 \iff 0\leq \theta \leq \pi.$$
Thus our new domain of integration is given by
$$\Gamma=\{(r,\theta,y)\in\mathbb{R}^3 : 0\leq r\leq 1 \land 0 \leq y \leq 2 \land 0\leq \theta \leq \pi\}.$$
We now turn to the actual integral. Using our change to cylindrical coordinates we get that
\begin{align*}
I&=\iiint_\Omega (x^2+y+z^2)^3 \mathrm{d}V =\iiint_\Gamma r(r^2\cos^2\theta+y+r^2\sin^2\theta)^3\mathrm{d}V \\
&=\iiint_\Gamma r(r^2(\cos^2\theta+\sin^2\theta)+y)^3\mathrm{d}V =\iiint_\Gamma r(r^2+y)^3\mathrm{d}V.
\end{align*}
Finally we can turn this triple integral into an iterated integral, yielding that
\begin{align*}
I &= \int_0^1 \int_0^2 \int_0^\pi r(r^2+y)^3 \mathrm{d}\theta~\mathrm{d}y~\mathrm{d}r = \pi\int_0^1 \int_0^2 r(r^2+y)^3 \mathrm{d}y~\mathrm{d}r \\
&= \pi\int_0^1 \biggl[ \frac{r(r^2+y)^4}{4} \biggr]_{y=0}^{y=2}~\mathrm{d}r = \pi\int_0^1 \frac{1}{4}\left(r(r^2+2)^4-r^9 \right)~\mathrm{d}r \\
&= \frac{\pi}{4}\int_0^1 \left(r(r^8+8r^6+24r^4+32r^2+16)-r^9 \right)~\mathrm{d}r \\
&= \frac{\pi}{4}\int_0^1 \left(r^9+8r^7+24r^5+32r^3+16r-r^9 \right)~\mathrm{d}r \\
&= \frac{\pi}{4}\int_0^1 \left(8r^7+24r^5+32r^3+16r \right)~\mathrm{d}r \\
&= \frac{\pi}{4} \biggl[ r^8+4r^6+8r^4+8r^2 \biggr]_0^1 \\
&= \frac{\pi}{4} (1+4+8+8) = \frac{21\pi}{4}.
\end{align*}
Thus, if I understood your region correctly, we have that $I=\frac{21\pi}{4}$