If $G$ is a domain (connected open set) that is symmetric with respect to the real axis, then $G\cap \mathbb{H} = \{ z\in G : \operatorname{Im} z > 0\}$ is also connected. So if $f\colon G_1 \to G_2$ is a holomorphic map between two domains symmetric with respect to the real axis, then $f(G_1\cap\mathbb{H})$ does not lie entirely above or below the real axis if and only if $f$ attains real values in $G_1\cap \mathbb{H}$. Thus, for the case of a biholomorphic $f$, it is sufficient to have $f(G_1\cap \mathbb{R}) = G_2 \cap \mathbb{R}$ to conclude that $f(G_1\cap \mathbb{H})$ is either the part of $G_2$ in the upper half-plane, or the part of $G_2$ in the lower half-plane.
If there is an $a\in G_1 \cap\mathbb{R}$ with $f(a) \in G_2 \cap \mathbb{R}$ and $f'(a) \in \mathbb{R}$, then $g(z) = \overline{f\left(\overline{z}\right)}$ is also a biholomorphic map $G_1 \to G_2$, with $g(a) = f(a)$ and $g'(a) = f'(a)$. Thus $h = g^{-1}\circ f$ is an automorphism of $G_1$ with fixed point $a$ and $h'(a) = 1$.
We obtain the conclusion if we can show that $h = \operatorname{id}_{G_1}$. For then $g = f$ shows that $f(G_1\cap\mathbb{R}) \subset G_2 \cap \mathbb{R}$, and the same argument for $f^{-1}$ shows the equality.
For that, we don't need any symmetry or other particular properties:
Theorem: Let $G\subset \mathbb{C}$ a domain, and $a\in G$. If $h$ is an automorphism of $G$ with $h(a) = a$ and $h'(a) = 1$, then $h = \operatorname{id}_G$.
Proof: If $G = \mathbb{C}$, then the automorphisms of $G$ are precisely the polynomials of degree $1$, and $h'(a) = 1$ then implies that $h$ is a translation $z \mapsto z+b$, and since $h(a) = a$, it follows that $b = 0$.
If $G = \mathbb{C}\setminus \{p\}$, then $h$ extends to an automorphism $\tilde{h}$ of the Riemann sphere $\widehat{\mathbb{C}}$, and either $\tilde{h}(p) = p,\, \tilde{h}(\infty) =\infty$, in which case $h = \operatorname{id}_{G}$ follows as above, or $\tilde{h}(p) = \infty$ and $\tilde{h}(\infty) = p$, in which case we have
$$\tilde{h}(z) = \frac{c}{z-p}+p$$
for some $c\in\mathbb{C}$. Then $h(a) = a$ yields $c = (a-p)^2$, and we obtain
$$h'(a) = -\frac{(a-p)^2}{(a-p)^2} = -1 \neq 1,$$
so that cannot happen: $h(a) = a$ and $h'(a) = 1$ implies $h = \operatorname{id}_{G}$ also for $G = \mathbb{C}\setminus \{p\}$.
If $\mathbb{C}\setminus G$ contains at least two points, then the family $\operatorname{Hol}(U,G)$ of holomorphic functions on any domain $U$ with values in $G$ is a normal family by Montel's big theorem. In particular - picking $U = G$ - the sequence $\left(h^n\right)_{n\in\mathbb{N}}$ of iterates of $h$ contains a locally uniformly convergent subsequence $\left(h^{n_k}\right)_{k\in\mathbb{N}}$. By Weierstraß' theorem, also the sequences $\left(\frac{d^m}{dz^m}h^{n_k}\right)_{k\in\mathbb{N}}$ of $m$-th derivatives converges locally uniformly for every $m\in\mathbb{N}$. If we already know that $h^{(r)}(a) = 0$ for $2 \leqslant r < m$ (which we vacuously do for $m = 2$), the Taylor expansion of $h$ about $a$ is
$$h(z) = a + (z-a) + c_m (z-a)^m + O\left((z-a)^{m+1}\right).$$
From that we find that the $n$-fold iterate of $h$ has the Taylor expansion
$$h^n(z) = a + (z-a) + n\cdot c_m (z-a)^m + O\left((z-a)^{m+1}\right)$$
by induction:
$$\begin{align}
h^{n+1}(z) &= h\left(a+(h^n(z)-a)\right)\\
&= a + \left(h^n(z)-a\right) + c_m \left(h^n(z)-a\right)^m + O\left(\left(h^n(z)-a\right)^{m+1}\right)\\
&= a + (z-a) + n\cdot c_m(z-a)^m + O\left((z-a)^{m+1}\right)\\
&\qquad + c_m \left((z-a) + O\left((z-a)^m\right)\right)^m + O\left((z-a)^{m+1}\right)\\
&= a + (z-a) + (n+1)c_m(z-a)^m + O\left((z-a)^{m+1}\right).
\end{align}$$
So
$$\frac{d^m}{dz^m}\Biggl\lvert_a h^n = n\cdot \frac{d^m}{dz^m}\Biggl\lvert_a h,$$
and that converges only if $h^{(m)}(a) = 0$. By induction on $m$, it follows that all derivatives of $h$ of order greater than one vanish in $a$, and hence $h = \operatorname{id}_{G}$.
In the plane, solving the Laplace equation with Dirichlet conditions is to find a function $u:\overline\Omega\subset\Bbb R^2\to\Bbb R$ such that: $$\begin{cases}\partial_{xx}u +\partial_{yy}u = 0, &(x,y)\in\Omega \\ u(x,y) = d(x,y), &(x,y)\in\partial\Omega \end{cases}$$
Verbally, we want a harmonic function $u$ defined in the domain $\Omega$, whose values coincide with those of some $\mathcal C^0 $ function $d$ on the border $\partial\Omega$. An easy version of the problem is when $\Omega = B_1(0)$ (the unit ball). Indeed, we are able to procure an analytic (in the sense that it's given by an explicit formula) solution, whose expression isn't important here, we just care that it exists.
Now say we wish to examine existence of solutions to the Laplace equation, where $\Omega$ is any ($\mathcal C^0$) simply connected region in $\Bbb R^2$ (except $\Bbb R^2$). You can probably see where this is going. Let $$\varphi:\Omega\to B_1(0)$$
be holomorphic (in the sense that the component functions $f,g$ make $f+ig$ holomorphic) with holomorhpic inverse (idem). Also, I'll require that $\varphi$ extend to $\partial\Omega$ continuously, and such that $\varphi(\partial\Omega) = \partial B_1(0)$. This is assured by a stronger version of the Riemann mapping theorem. Now consider the Dirichlet problem $$\begin{cases}\partial_{xx}u +\partial_{yy}u = 0, &(x,y)\in B_1(0) \\ u(x,y) = (d\circ \varphi^{-1})(x,y), &(x,y)\in\partial B_1(0) \end{cases}$$
This is exactly of the type that I first mention, so it has a solution $v$. I now claim that $v\circ \varphi$ is a solution to the Dirichlet problem $$\begin{cases}\partial_{xx}u +\partial_{yy}u = 0, &(x,y)\in \Omega \\ u(x,y) = d(x,y), &(x,y)\in\partial\Omega \end{cases}$$
The boundary condition is immediate, since on $\partial B_1(0)$, $v=d\circ \varphi^{-1}$. Therefore on $\partial\Omega$, $v\circ\varphi = d\circ\varphi^{-1}\circ\varphi = d$.
For the PDE, we use that since $v$ is harmonic (i.e. it satisfies the Laplace equation), there exists a harmonic conjugate $w$ such that $v+iw$ is holomorphic. Therefore $(v+iw)\circ\varphi$ (now interpreting $\varphi$ as complex) is holomorphic, and therefore its real part $v\circ\varphi$ is harmonic*, i.e. it satisfies $(\partial_{xx}+\partial_{yy})(v\circ\varphi) = 0$.
By the Riemann mapping theorem, we were able to solve the Laplace equation with Dirichlet conditions over a much wider range of domains, using only that we know a solution in the unit ball. Of course there are much more potent ways of approaching this particular equation, that give even more general results, but I think it's a nice example of an application regardless.
*It's easy to check that if $f+ig$ is holomorphic, then $f,g$ are harmonic, by using the Cauchy-Riemann equations
Best Answer
I've found following resource:
http://www.dancalloway.com/lets_have_a_word/riemann-mapping-theorem-explained
It's really very good and fitted my needs well.