This theorem seems almost magical. The algebraic derivation doesn't really provide any insight into why it works.
So could someone give me a geometric interpretation of it?
This: Geometrical Interpretation of the Cauchy-Goursat Theorem? is a similar question but uses another, incorrect, definition of the theorem.
On
The Cauchy(-Goursat) theorem deals with arbitrary analytic functions. Your formula in red is a simple particular instance.
When an analytic function $f:\>\Omega\to{\mathbb C}$ has a primitive $F:\>\Omega\to{\mathbb C}$ then automatically the integral $\int_\gamma f(z)\>dz$ is zero for all closed curves $\gamma \subset \Omega$, because the integral of $f$ along any arc from $p$ to $q$ is equal to the difference $F(q)-F(p)$.
Since $f(z):=(z-z_0)^n$, $\>n\ne-1$, has the primitive $F(z)={1\over n+1}(z-z_0)^{n+1}$ in $\Omega:={\mathbb C}\setminus\{z_0\}$ this principle takes care of the claim for all $n\ne-1$.
When dealing with the remaining case I shall assume $z_0=0$; so it's about the integral $$\int_C{dz\over z}\ .$$ By assumption the curve $C$ (assumed smooth for simplicity) goes once around the origin. Therefore it has a parametric representation of the form $$C:\quad t\mapsto z(t):=r(t)e^{i\phi(t)}\qquad(0\leq t\leq1)$$ with $r(t)>0$ for all $t$, $r(1)=r(0)$, and $$\phi(1)-\phi(0)=2\pi\ .$$ Note that $\phi(\cdot)$ need not be monotone. It follows that $$\int_C{dz\over z}=\int_0^1{z'(t)\over z(t)}\>dt=\int_0^1\left({r'(t)\over r(t)}+i\phi'(t)\right)\>dt=\left(\log\bigl(r(t)\bigr)+i\>\phi(t)\right)\biggr|_0^1=2\pi\>i\ .$$
Let's take the path integral around the origin at radius 1, for convenience. On this circle, we have $1/z = \bar{z}$, the complex conjugate of $z$. If $f(z) = 1/z$ then $z \cdot f(z) = 1$ is real, and if $w$ is any complex number on the same ray from the origin as $z$, then $w \cdot f(z)$ is also real. So in particular any infinitesimal step along the radius $z \to z(1 + \mathrm{d}r)$, with $\mathrm{d}r$ real, would make a purely real contribution to an integral of $f(z)\mathrm{d}z$.
Similarly (because $f$ is a conformal map) any infinitesimal step around the circumference of of the unit circle makes a purely imaginary contribution to the integral. As z goes about the circle anticlockwise, $f(z) = \bar{z}$ goes clockwise at just the right speed to keep up and make the contribution constant. These imaginary parts add up over the length of the circumference to $2 \pi i$. (Note that the scaling at this radius is exactly 1).
Now let's take $f(z) = 1/z^n$ for some higher $n$. Now the argument of $f(z)$ rotates at a higher frequency as we go around the unit circle; instead of all the contributions being (positive) imaginary and adding up, they point evenly in all directions. Indeed we get $n-1$ contributions in each direction from zero. So all the contributions cancel out and the integral is zero.