The motivation is that you're trading the sum of all finite residues for a single "residue at infinity." For example, if you had only a finite number of poles in the plane for a holomorphic function $F$, and $\Gamma$ enclosing those poles in its interior, then the sum of the residues would be
$$
\frac{1}{2\pi i}\oint_{\Gamma}F(\lambda)d\lambda = -\frac{1}{2\pi}\oint_{1/\Gamma}F(1/\mu)\frac{1}{\mu^2}d\mu.
$$
The negative would cancel the new negative orientation of $1/\Gamma$, and the end result of the integration as you let $\Gamma$ expand without bound in, say, a circle, would be
$$
\lim_{\mu\rightarrow 0}F(1/\mu)(1/\mu)= \lim_{\lambda\rightarrow \infty}\lambda F(\lambda).
$$
There are lots of issues in how that limit is achieved, but that's the basic idea. This idea works very nicely for matrices and operators, too. Fredholm was the pioneer in this type of analysis, and his work fueled earliest forms of Spectral Theory through the use of the resolvent operator $R(\lambda)=(L-\lambda I)^{-1}$. (Resolvent was a term coined by Fredholm, who also was the first to define a linear operator.) For example, if you have an $N\times N$ selfadjoint matrix $L$, then you can show that $(\lambda I-L)^{-1}$ has simple poles at the eigenvalues, and the residue at such a pole is the projection onto the eigenspace associated with that eigenvalue. Then completeness of the eigenfunctions is a consequence of the fact that the sum of these residues is
$$
\lim_{\lambda\rightarrow\infty}\lambda (\lambda I-L)^{-1} = \lim_{\lambda\rightarrow\infty}\frac{\lambda}{\lambda I-L}=I.
$$
And you really can make this rigorous because you can show the above limit is $I$ for any $N\times N$ matrix. Furthermore, regardless of the matrix, you can show that the residue at an eigenvalue is a projection onto the space spanned by the Jordan blocks associated with the eigenvalue. So you have completeness due to an unusual conservation law associated with holomorphic functions. By looking on the Riemann sphere, you see that you can trade the finite residues for the residue at $\infty$. For normal and selfadjoint matrices, all poles are first order the and the residues are the projections onto the corresponding eigenspaces. And the Complex Analysis trick then shows that the sum of all of the projections onto eigenvalues of a normal matrix must sum to the identity $I$, thereby proving that the eigenvectors form a basis.
This trick also works for more general operators, such as the differentiation operator $Lf = \frac{1}{i}\frac{d}{dx}$ on $L^2[-\pi,\pi]$ where the domain is chosen to consist of continuously differentiable periodic functions with $f' \in L^2$ (absolutely continuous is even better.) The resolvent can be computed directly for this operator by solving $(L-\lambda I)g=f$ for $g$ as a function of $\lambda$. This is done by assuming a given $f$, and solving a first order ODE for $g$ as a function of $\lambda$:
$$
\frac{1}{i}g'-\lambda g = f \\
g(0)=g(2\pi)
$$
It is my understanding that Cauchy came up with the proof given by Whittaker by considering this equation. To solve this equation multiply by $i$ and then by an integrating factor $e^{-i\lambda t}$ to obtain
$$
e^{-i\lambda t}g'-i\lambda e^{-i\lambda t}g = ie^{-i\lambda t}f \\
\frac{d}{dt}(e^{-i\lambda t}g)=ie^{-i\lambda t}f \\
e^{-i\lambda t}g(t) = i\int_{0}^{t}e^{-i\lambda s}f(s)ds + C \\
g(t) = ie^{i\lambda t}\int_{0}^{t}e^{-i\lambda s}f(s)ds+Ce^{i\lambda t}
$$
The constant $C$ is determined by requiring periodicity:
$$
C = g(0) = g(2\pi)=ie^{2\pi i\lambda}\int_{0}^{2\pi}e^{-i\lambda s}f(s)ds+Ce^{2\pi i\lambda} \\
C(1-e^{2\pi i\lambda})=ie^{2\pi i\lambda}\int_{0}^{2\pi}e^{-i\lambda s}f(s)ds \\
C = \frac{ie^{2\pi i\lambda}}{1-e^{2\pi i\lambda}}\int_{0}^{2\pi}e^{-i\lambda s}f(s)ds
$$
Therefore, $g(t,\lambda)$ is given by
$$
g(t,\lambda)=\left(i\int_{0}^{t}f(s)e^{-i\lambda s}ds
+\frac{ie^{2\pi i\lambda}}{1-e^{2\pi i\lambda}}\int_{0}^{2\pi}e^{-i\lambda s}f(s)ds\right) e^{i\lambda t}
$$
Notice that the residues of this expression are negatives of the projections onto the one-dimensional eigenspaces associated with $e^{int}$.
$$
R_n = -\frac{1}{2\pi}\int_{0}^{2\pi}f(s)e^{-ins}ds e^{int}
$$
(The negative is because of using $(L-\lambda I)^{-1}$ instead of $(\lambda I-L)^{-1}$.)
You can write this in a more symmetric form simply by splitting the integral over $[0,2\pi]$ into integrals over $[0,t]$ and $[t,2\pi]$, and that's how the analysis is carried out for $\lim_{\lambda\rightarrow\infty}\lambda g(t,\lambda)$, after noting that the $g$ associated with a constant $f \equiv C$ is easily found to be $-C/\lambda$ because a constant function is periodic, and $(\frac{1}{i}\frac{d}{dt}-\lambda)\frac{-C}{\lambda}=C$. (I believe Whittaker's analysis stems from the use of the symmetric form.)
Reference: E. C. Titchmarsh, Eigenfunction Expansions Associated with Second Order Differential Equations
Titchmarsh proves a general expansion theorem that includes the ordinary Fourier case in the first 20 pages of the first chapter. Titchmarsh was student of G. H. Hardy, and pioneered much of the rigorous pointwise analysis for this subject.
Further Detail: Rewrite $g(x,\lambda)=R(\lambda)f$ as
\begin{align}
R(\lambda)f & = \frac{i}{1-e^{2\pi i\lambda}}\left\{
\int_{0}^{x} e^{i\lambda (x-t)}f(t)\,dt -
\int_{x}^{2\pi}e^{i\lambda(2\pi-(t-x))}f(t)\,dt\right\} \\
& = \frac{i}{1-e^{-2\pi i\lambda}} \left\{
-\int_{0}^{x} e^{-i\lambda(2\pi-(x-t))}f(t)\,dt
+ \int_{x}^{2\pi}e^{-i\lambda(t-x)}f(t)\,dt \right\}
\end{align}
The first form is convenient for examining the resolvent for $\Im\lambda > 0$, and the second is convenient for $\Im\lambda < 0$. When examining the integral of the resolvent on a circle of half-integer radius $|\lambda|=N+1/2$, the function $1/(1-e^{2\pi i\lambda})$ is uniformly bounded by a constant $M$ for $\Im\lambda \ge 0$, and the exponentials in the integrals are well-behaved. A similar analysis may be carried out using the second form for $\Im\lambda \le 0$.
More to come later ...
The usual proof of Stokes's theorem requires the differential form to be at least of class $C^1$ (i.e., differentiable everywhere with continuous first partial derivatives). There are versions that work with weaker hypotheses, but just assuming the existence of derivatives is definitely not enough.
Your proof of Cauchy's formula works fine as long as you assume that $f$ is of class $C^1$ and satisfies the Cauchy-Riemann equations. But the Cauchy-Goursat argument is needed if you only assume that $f'(z)$ exists for all $z\in\Omega$, because then Stokes’s theorem doesn't apply.
Best Answer
Let $X$ be a compact Riemann surface and $\omega$ a smooth one-form on $X-p_1,\ldots,p_N$, let $U$ be an open $\subset X$ such that : $U$ contains all the $p_n$, it is homeomorphic to a simply connected domain in $\Bbb{C}$ and $U\cup\partial U=X$.
From here see $\partial U$ not as a subset of $X$ but as a closed-loop in $X$.
Topologically $X$ is obtaining by gluing different pieces of $\partial U$ together, that is $\partial U = \bigcup_{j=1}^J \gamma_j$ and each $\gamma_j$ is glued with some $\gamma_{\sigma(j)}$.
The complex topology induces an orientation, and as curves in $X$,
$\gamma_{\sigma(j)}$ must be $\gamma_j$ traversed in opposite direction because otherwise the points "on the left side of $\gamma_j$" would appear twice in $U$.
Thus $\partial U$ is homologous to $0$
$$2\int_{\partial U} \omega=\sum_j \int_{\gamma_j}+\int_{\gamma_{\sigma(j)}} \omega=\sum_j 0=0$$
And hence if $\omega$ is a meromorphic one-form ($\omega=gdf$ for two meromorphic functions) then
$$2i\pi\sum Res_{p_n}(\omega) = \int_{\partial U} \omega=0$$