Why are the integral curves of this vector field easy to find?

Indentifying $(x,y)$ with $z=x + y \, \mathrm{i}$ the differential equation reads $$\frac{\partial z}{\partial t}=z^2.$$ Rewriting this we find $$\frac{-\partial z^{-1}}{\partial t} = \frac{\partial z}{\partial t} z^{-2} = 1$$ and so $$-z^{-1}= t+z_0 \textrm{ or } z = -\frac1{t+z_0}$$ for some fixed initial value $z_0 \in \mathbb{C}$. For $t \in \mathbb{R}$ these describe circular trajectories as expected.

Alternative 1: transform the vector field under the inversion

$$\sigma(x, y) = \left(\frac{x}{x^2+y^2}, \frac{y}{x^2+y^2}\right).$$

This map is an involution ($\sigma^2 = \mathrm{id}$) and its Jacobian is $$D\sigma_{(x,y)}=\frac{1}{(x^2+y^2)^2}\begin{pmatrix} y^2-x^2 & -2xy \\ -2xy & x^2-y^2 \end{pmatrix}.$$

Therefore it maps the element $$\begin{pmatrix}x^2 - y^2 \\ 2xy \end{pmatrix}$$ at point $(x, y)$ to the element

$$\frac{1}{(x^2+y^2)^2}\begin{pmatrix} y^2-x^2 & -2xy \\ -2xy & x^2-y^2 \end{pmatrix} \begin{pmatrix}x^2 - y^2 \\ 2xy \end{pmatrix} = \begin{pmatrix} -1 \\ 0 \end{pmatrix}$$ at point $\sigma(x,y)$. The trajectories of the transformed vector field are easy to find. Then apply the inversion $\sigma$ to find the trajectories of the original vector field.

Alternative 2: Since you know that the trajectories are circular start with the ansatz $$\begin{cases} x = -R \sin \alpha \\ y = R + R \cos \alpha \end{cases}$$ for some function $\alpha$ of $t$. Then the given differential equations are consistent and both lead to $$\frac{\partial \alpha}{\partial t} = 2R \, (1 + \cos\alpha).$$ One solution to this equation is $\alpha(t) = 2 \arctan(2R \, t)$. Substituting this in the expressions for $x$ and $y$ gives the trajectory

$$\left(\frac{-t}{t^2 + (2R)^{-2}}, \frac{(2R)^{-1}}{t^2+(2R)^{-2}}\right).$$

Note that the system implies $$y'(x)=\frac{2xy(x)}{x^2-y(x)^2}.$$ (This is the approach mentioned by @Sameh Shenawy in the comments.) This can be solved by standard methods. For example, setting $y(x)=xf(x)$ we arrive at the separable differential equation $$f'(x) = \frac{1}{x} \frac{f(x)(1+f(x)^2)}{1-f(x)^2}.$$ The relevant integral may be done by a partial fraction expansion.

In my opinion the method put forward by @WimC is much preferable.

Addendum

The gory detail ... \begin{align*} f'(x) &= \frac{1}{x} \frac{f(x)(1+f(x)^2)}{1-f(x)^2} \\ \frac{1-f^2}{f(1+f^2)}df &= \frac{dx}{x} \\ \frac{(1+f^2)-2f^2}{f(1+f^2)} df &= \frac{dx}{x} \\ \int\left(\frac{1}{f} - \frac{2f}{1+f^2}\right)df &= \int\frac{dx}{x} \\ \ln f - \ln(1+f^2) &= \ln\frac{x}{2r} \\ \ln\frac{f}{1+f^2} &= \ln\frac{x}{2r} \\ \frac{f}{1+f^2} &= \frac{x}{2r} \\ \frac{y/x}{1+y^2/x^2} &= \frac{x}{2r} \\ \frac{y}{x^2+y^2} &= \frac{1}{2r} \\ x^2+y^2-2ry &= 0 \\ x^2 + y^2 - 2ry + r^2 &= r^2 \\ x^2 + (y-r)^2 &= r^2. \end{align*} (The integration constant was chosen to have the form $-\ln 2r$ for convenience.)