Amazingly, it seems that the answer is yes:
Mioduszewski, J.
On two-to-one continuous functions. (Russian summary)
Bull. Acad. Polon. Sci. Sér. Sci. Math. Astronom. Phys. 9 1961 129--132.
The author announces results concerning two-to-one functions $f$ on a locally compact separable space $X$, proofs of which appear in Rozprawy Mat. 24 (1962), 1--41. Let $\phi$ be the (discontinuous) involution defined by $\varphi(x)=f^{-1}f(x)-x$. A result of the reviewer [Duke Math. J. 10 (1943), 49--57; MR0008697 (5,47e)] asserts that if $X$ is a compact manifold or $f$ is closed and $X$ is a locally compact manifold, then the investigation of $\phi$ is equivalent to the investigation of a continuous involution. The author calls a point $x\in X$ pseudo-Euclidean if it has a neighborhood $H$ such that the closure of the component of $x$ in $H$ is a Euclidean solid sphere. The principal theorem asserts that if $x$ is a pseudo-Euclidean point with $K$ as the solid sphere of the definition, and if $\psi=\varphi|K$, that $\lim\text{}\sup_{y\rightarrow x}\psi(y)=x\bigcup\varphi(x)$ is impossible. This yields an extension of the result of the reviewer quoted above. The author indicates the existence of a plane simply connected domain $G$ whose boundary is an irreducible cut of the plane into two domains and such that there exists a two-to-one mapping defined on $\overline G$. This is in contrast to the result of Roberts [ibid. 6 (1940), 256--262; MR0001923 (1,319d)], which asserts the non-existence of two-to-one mappings defined on two-cells. The existence of two-to-one mappings defined on Euclidean spaces $E^n$, $n\geq 2$, is shown. However, the question of the existence of two-to-one mappings defined on $n$-cells, $n>3$, remains open. [MathSciNet review by P. Civin.]
I can't access this paper, so I can't say anything about the construction. It would be nice to see some corroboration for this result and/or a more (physically) accessible contemporary treatment.
Addendum: Petya's response gives a link to the paper, from which one can see that the function is essentially defined in terms of the involution $\iota$, so it is not immediately clear what the codomain is or whether it can be embedded in $\mathbb{R}^2$.