I think we could exploit the fact that the integrand is a gradient of the fundamental function for a Poisson equation.
Notice:
$$
\nabla \big(\ln((x+1)^2+y^2)\big) = 2\left(\frac{x+1}{(x+1)^2+y^2},
\frac{y}{(x+1)^2+y^2}\right)
$$
Denote the square as $\Omega$, the infinitesimal tangent of the curve is $(dx,dy)$, then rotate it clockwise by $\pi/2$ we get the normal, and the normal is $(dy,-dx)$. So the original integral $I$ becomes:
$$
I = -\frac{1}{2}\int_{\partial \Omega} \nabla \big(\ln((x+1)^2+y^2)\big) \cdot \mathbf{n} \,ds
$$
Moreover recall the fundamental solution of Poisson equation in 2D is $\displaystyle \phi = \frac{1}{2\pi}\ln|\mathbf{x}-\mathbf{x}_0|$ which solves:
$$
\Delta \phi = \delta(\mathbf{x}-\mathbf{x}_0)
$$
in $\mathbb{R}^2$. Let $\mathbf{x}_0 = (-1,0)$, then $u =\ln\big((x+1)^2+y^2\big) = 2\ln|\mathbf{x}-\mathbf{x}_0|$, and $u$ solves:
$$
\Delta u = 4\pi \delta(\mathbf{x}-\mathbf{x}_0)
$$
Now uses Divergence theorem:
$$
I = -\frac{1}{2}\int_{\partial \Omega} \nabla u \cdot \mathbf{n} \,ds = -\frac{1}{2} \int_{\Omega} \Delta u = -\frac{1}{2} \int_{\Omega} 4\pi \delta(\mathbf{x}-\mathbf{x}_0) = -2\pi
$$
Above method looks like a little bit cheating though.
EDIT: Guess this is a standard Calculus III problem, so we shouldn't use either Dirac delta nor Cauchy integral formula. Thinking of the Cauchy integral formula reminded me a trick we used in Calculus III (maybe it was too long and slipped from my mind, now I recollected how we did this kind of problem):
When the enclosed $\Omega$ by the curve $C$ has a singularity, i.e., $\displaystyle \frac{\partial Q}{\partial x}$ or $\displaystyle \frac{\partial P}{\partial y}$ is not continuous, we could just cut a hole with a small radius $r$ centered at this singularity , in this example it is $(-1,0)$, such that outside this small disk of radius $r$ there is no singularity (See the picture below).

Denote $C'$ the boundary of this disk $B$, rotating clockwisely. Then:
$$
\oint_{C} \frac{ydx-(x+1)dy}{x^2+y^2+2x+1} + \oint_{C'} \frac{ydx-(x+1)dy}{x^2+y^2+2x+1} =
\int_{\Omega\backslash B}\left\{ -\frac{\partial}{\partial x}\left(\frac{x+1}{(x+1)^2+y^2}\right) - \frac{\partial}{\partial y}\left(\frac{y}{(x+1)^2+y^2}\right)\right\} dx dy = 0
$$
Therefore:
$$
\oint_{C} \frac{ydx-(x+1)dy}{x^2+y^2+2x+1} =- \oint_{C'} \frac{ydx-(x+1)dy}{x^2+y^2+2x+1}
$$
Now we could just parametrize the curve $C'$ using $t$, by letting $x+1 = r\cos t$, $y = r\sin t$, so $dx = -r\sin t\, dt$, $dy = r \cos t \,dt$. The integral becomes:
$$
\oint_{C'}\frac{y }{(x+1)^2+y^2}dx-\frac{(x+1)}{(x+1)^2+y^2}dy
= \int^{2\pi}_0 \left(\frac{r\sin t\, }{r^2\cos^2 t + r^2\sin^2 t} (-r\sin t)
- \frac{r\cos t\,}{r^2\cos^2 t + r^2\sin^2 t}r \cos t \right)dt = -2\pi
$$