The best way to understand how or why the author made that specific choice for $\delta$ and $M$ a priori, is to commence a proof using only the information at hand and see what develops. As you will see, other choices are possible and need not be determined at the start by some flash of insight.
The objective is to find $\delta > 0$ such that for all $(x,y) \in W$ we have $|xy - x_0y_0| < \epsilon.$
You have to anticipate that the determination of a suitable $\delta$ will require some estimation and bounding. The tools we have immediately at our disposal are $$(x,y) \in [x_0-\delta,x_0 + \delta] \times [y_0-\delta,y_0 + \delta],$$ which implies
$$|x - x_0| \leqslant \delta, \quad |y - y_0| \leqslant \delta,$$
along with,
$$\quad |x| \leqslant |x_0| + |x - x_0| \leqslant |x_0| + \delta, \\ \,\,\,\,|y| \leqslant |y_0| + |y - y_0| \leqslant \,\,|y_0| + \,\delta $$
The last two lines follow from the reverse triangle inequality or, if you prefer, the usual triangle inequality with $|x| = |x_0 + (x - x_0) | \leqslant |x_0| + |x-x_0|$.
Since we have only simple inequalities involving $x$ and $y$ separately, we make the decomposition $|xy- x_0y_0| \leqslant |x - x_0||y| + |x_0||y - y_0|$ and then obtain using those inequalities,
$$\tag{*}|xy- x_0y_0| \leqslant |x - x_0||y| + |x_0||y - y_0| \leqslant \delta(|y_0| + \delta) + \delta|x_0|\\ \leqslant 2\delta \max(|x_0|,|y_0| + \delta)$$
We want to find $\delta$ such that the RHS of (*) is smaller than $\epsilon$, but at this point the procedure will be messy. When $|y_0| > |x_0|$ we are left with a quadratic equation, and to examine the other case where $|x_0| > |y_0| + \delta$ we are faced with an iterative process.
However, if we are not concerned with finding the largest possible $\delta$, then we can set a constraint $\delta < \alpha$ where $\alpha > 0$ is arbitrary and replace $\max(|x_0|, |y_0| + \delta)$ with $M = \max(|x_0|, |y_0| + \alpha)$, to obtain
$$|xy- x_0y_0| \leqslant 2\delta M$$
If both $\delta \leqslant \alpha$ and $\delta \leqslant \epsilon/(2M)$, that is $\delta \leqslant \min(\alpha,\epsilon/(2M))$, then $|xy-x_0y_0| \leqslant \epsilon$.
The author happened to like another choice, $M = \max(|x_0|+1,|y_0| + 1)$, perhaps because it looks "cleaner", and knew this would work in advance at the start of the proof.
Finally, you are correct in that the statement should be $xy$ belongs to a neighborhood of $x_0y_0$ of radius $\epsilon$. Here neighborhood means closed ball with center $x_0y_0$ and radius $\epsilon$.
Best Answer
Choose $\delta = \varepsilon/2$. Then, for all $(x,y)$ such that $\| (x,y)-(x_0,y_0) \| < \delta$ we have that $$\max \{ |x-x_0|,|y-y_0| \} \leq \| (x-x_0,y-y_0) \| = \| (x,y)-(x_0,y_0) \| < \frac{\varepsilon}{2}$$ which means $|x-x_0|<\varepsilon/2$ and $|y-y_0|<\varepsilon/2$. Thus $$|f(x,y)-f(x_0,y_0)|=|(x+y)-(x_0+y_0)|\leq |x-x_0|+|y-y_0|<\varepsilon.$$