Here is an explanation based on the Euler-Maclaurin summation formula.
(Or rather, since we'll only ever need two terms of the Euler-Maclaurin summation, it's really more or less just "the trapezoid rule".)
I think it's a good explanation because it sticks to the general structure of the argument outlined in the question, and its "18th-century-friendly" spirit.

First, let's review how to apply Euler-Maclaurin to $\zeta$.
We have
$$
\begin{align}
\zeta(s) &=
\sum_{n=1}^N \frac{1}{n^s}
+ \sum_{n=N+1}^{\infty} \frac{1}{n^s} \\
&=
\sum_{n=1}^N \frac{1}{n^s}
+ \int_N^{\infty} \frac{1}{x^s} \, dx
- \frac12 \frac{1}{N^s}
+ \mathrm{Error} \\
&=
\sum_{n=1}^N \frac{1}{n^s}
+ \frac{1}{s-1} \frac{1}{N^{s-1}}
- \frac12 \frac{1}{N^s}
+ \mathrm{Error}
\tag{1}. \\
\end{align}
$$
We can say more about the error term later - for now, all we need to know is that $\mathrm{Error} = O ( \frac{1}{N^{\operatorname{Re}(s)+1}} )$ as $N \to \infty$, for each $s$.
The key thing to note is that, even though the computation initially required $\operatorname{Re}(s)>1$ in order to be valid, in fact the "equation"
$$
\zeta(s) =
\sum_{n=1}^N \frac{1}{n^s}
+ \frac{1}{s-1} \frac{1}{N^{s-1}}
- \frac12 \frac{1}{N^s}
+ O \left( \frac{1}{N^{\operatorname{Re}(s)+1}} \right)
\tag{2}
$$
has a unique constant solution $\zeta(s)$, for each $s$ with $\operatorname{Re}(s) > -1$ and $s \ne 1$.
This is one way to *define* $\zeta(s)$ for all such $s$.
In particular, we can plug in $0$ and immediately find that $\zeta(0) = -\frac12$.

What about the function $f(z) = -\frac{\pi z}{2} \cot(\pi z)$? We can also apply Euler-Maclaurin to $f$ in the same way:
$$
\begin{align}
f(z) &=
-\frac12
+ \sum_{n=1}^N \frac{z^2}{n^2-z^2}
+ \sum_{n=N+1}^{\infty} \frac{z^2}{n^2-z^2} \\
&=
-\frac12
+ \sum_{n=1}^N \frac{z^2}{n^2-z^2}
+ \int_N^{\infty} \frac{z^2}{x^2-z^2} \, dx
- \frac12 \frac{z^2}{N^2-z^2}
+ O \left( \frac{1}{N^3} \right)
\tag{3} \\
\end{align}
$$
as $N \to \infty$, for each $z$.
The right-hand side of (3), without the error term, is what we'll call $f_N(z)$; we can simplify it to
$$
f_N(z) =
\sum_{n=1}^N \frac{z^2}{n^2-z^2}
+ z \cdot \frac12 \left(
\log\left(1+\frac{z}{N}\right)
- \log\left(1-\frac{z}{N}\right)
\right)
- \frac12 \frac{N^2}{N^2-z^2}
\tag{4}.
$$
I should emphasize that (4) is just a somewhat more elaborate variant of $-\frac12 + \sum_{n=1}^{\infty} \frac{z^2}{n^2-z^2}$, much like how (2) is just a somewhat more elaborate variant of $\sum_{n=1}^{\infty} \frac{1}{n^s}$.

Now when we expand (4) into power series, we recognize the expression from (2) as the coefficients, and we recognize that we can make the sum run from $k=0$ to $\infty$, not just $k=1$ to $\infty$:
$$
\begin{align}
f_N(z) &=
\sum_{n=1}^N \sum_{k=1}^{\infty} \frac{z^{2k}}{n^{2k}}
+ \sum_{k=1}^{\infty} \frac{1}{2k-1} \frac{1}{N^{2k-1}} z^{2k}
- \frac12 \sum_{k=0}^{\infty} \frac{z^{2k}}{N^{2k}} \\
&=
\sum_{k=0}^{\infty} \left(
\sum_{n=1}^N \frac{1}{n^{2k}}
+ \frac{1}{2k-1} \frac{1}{N^{2k-1}}
- \frac12 \frac{1}{N^{2k}}
\right) z^{2k}.
\tag{5} \\
\end{align}
$$
This is exactly what we want: take the limit as $N \to \infty$ to conclude that $f(z) = \sum_{k=0}^{\infty} \zeta(2k) z^{2k}$.

(To be rigorous in the final step, we'd need to be more precise about the error term in (1). The actual bound we get from Euler-Maclaurin is
$
\lvert \mathrm{Error} \rvert
\le \mathrm{constant} \cdot \int_N^{\infty} \left\lvert \frac{s(s+1)}{x^{s+2}} \right\rvert \, dx
$
for all $N$ and all $s$ such that $\operatorname{Re}(s) > -1$ and $s \ne 1$.
This lets us control the size of the difference $\sum_{k=0}^{\infty} \zeta(2k) z^{2k} - f_N(z)$.)

This proves that all the even-power coefficients of the power series of $-\frac{\pi z}{2} \cot(\pi z)$ are the corresponding values of $\zeta$, including $\zeta(0)$ as the constant term.

Divergent series, discusses this, I think in the last chapter. Another place I think I saw this is P. Cartier's contribution in the volumeFrom Number Theory to Physics, Springer 1992. (I don't have any of these sources in front of me.) $\endgroup$ – Liviu Nicolaescu Nov 29 '14 at 12:39