We need to investigate if there exists constant $\lambda \in [1,1)$ such that
$$
\forall p_1, p_2 \in \mathbb{R}^2 \: : \: \|\Phi(p_1) - \Phi(p_2)\|_2 \leq \lambda \|p_1-p_2\|_2$$
and a constant $\mu \in [0,1)$
$$
\forall p_1, p_2 \in \mathbb{R}^2 \: : \: \|\Phi(p_1) - \Phi(p_2)\|_\infty \leq \mu \|p_1-p_2\|_\infty$$
Let therefore $p_1 = (x_1,y_1)$ and $p_2 = (x_2, y_2)$ be two points in $\mathbb{R}^2$. We must control the difference between $\Phi(p_1)$ and $\Phi(p_2)$. When $n=1$ we can typically apply the mean value theorem. In higher dimensions a rewrite it is required:
Define the auxiliary function $g : [0,1] \rightarrow \mathbb{R}^2$ by
$$ g(t) = \Phi(tp_2+(1-t)p_1).$$
By the chain rule, this function is differentiable and the derivative is
$$ g'(t) = D\Phi(tp_2+(1-t)p_1)\cdot(p_2-p_1) dt.$$
Here $D\Phi$ denotes the Jacobian of $\Phi$ and $\cdot$ is the Euclidian inner product. Moreover $g(0)=\Phi(p_1)$ and $g(1) = \Phi(p_2)$. It follows that
$$ \Phi(p_2) - \Phi(p_1) = g(1) - g(0) = \int_0^1 g'(t) dt = \int_{0}^1 D\Phi(tp_2+(1-t)p_1)\cdot(p_2-p_1) dt.$$
This implies
$$ \| \Phi(p_2) - \Phi(p_1) \| \leq \int_{0}^1 \| D\Phi(tp_2+(1-t)p_1) \| \cdot \|(p_2-p_1) \|dt. $$
regardless of the choice of norm. We conclude that it is worthwhile to investigate the size of
$$ C_2 = \int_{0}^1 \| D\Phi(tp_2+(1-t)p_1) \|_2 dt$$
as well as
$$ C_\infty = \int_{0}^1 \| D\Phi(tp_2+(1-t)p_1) \|_\infty dt$$
Now, if the Jacobians could be bounded by a number strictly less than unity, we would be done. With $p = (x,y)$ we have
$$ D\Phi(p) = \frac{1}{2} \begin{bmatrix} \frac{1}{4} \cos(x) & 1 \\ 1 & \cos(y) \end{bmatrix}. $$
The Jacobian is real and symmetric, so the 2-norm is determined by the eigenvalues. By Gershgorin's circle theorem, it follows that
$$ \|D\Phi(p)\|_2 \leq 1.$$
By direct definition of the infinity norm we have
$$ \|D\Phi(p)\|_\infty \leq 1.$$
This certainly does not preclude the possibility that $\Phi$ is a contractive mapping, but more work is required.
To that end, we return to the expression for $D\Phi(p)$. For simplicity we consider $2 D\Phi(p)$. The characteristic polynomial of $2 D\Phi(p)$ is
$$ p(\lambda) = \lambda^2 - \left(\frac{1}{4} \cos(x) + \cos(y) \right)\lambda + \left(\frac{1}{4}\cos(x)\cos(y)-1\right)$$
The pressing issue is bound the eigenvalues away from $\pm 2$. By the triangle inequality a root $\lambda$ must necessarily satisfy
$$ |\lambda|^2 \leq \frac{5}{4} |\lambda| + \frac{5}{4}.$$
This puts a hard limit on the size of $|\lambda|$, specifically
$$|\lambda| \leq \frac{5 + \sqrt{105}}{4} < 2.$$
We can in fact conclude that $\Phi$ is a contractive mapping with respect to the $2$-norm.
Currently, I have no good ideas about the case of the infinity norm.
UPDATE: Let us explore what can be done by direct calculation with $\Phi$. We have
$$ \Phi(p_2) - \Phi(p_1) = \frac{1}{2} \begin{bmatrix} \frac{1}{4} \left(\sin(x_2) - \sin(x_1)\right) + y_2 - y_1 \\ x_2-x_1 + \sin(y_2) - \sin(y_1) \end{bmatrix} $$
The following inequality should be understod in the componentwise sense:
$$|\Phi(p_2) - \Phi(p_1)| \leq \frac{1}{2} \begin{bmatrix} \frac{1}{4}|x_2 - x_1| + |y_2-y_1| \\ |x_2 - x_1| + |y_2 - y_1| \end{bmatrix}$$
It follows
$$ \|\Phi(p_2) - \Phi(p_1)\|_\infty \leq \|p_2-p_1\|_\infty $$
Moreover, exploring the choice of $x_1=y_1=h$ and $x_2=y_2=2h$ and allowing $h$ to tend to zero will reveal that this is the best possible bound. We conclude that $\Phi$ is not a contraction with respect to the infinity norm.
