Any linear fractional transformation transforming the real axis to itself can be written in terms of reals?

Given $\frac{az+b}{cz+d}$, you can trivially assume that $a$ is real. Either $a=0$, or you can multiply all four numbers by $\overline{a}$. (or divide all four by $a$)

Now that $a$ is real, $c$ can be assumed real too. If $a$ were zero, apply the same trick to $c$. ($c$ is not zero, or else the transformation is noninvertible. So multiply all coefficients by $\overline{c}$ or divide them all by $c$.) If $a$ were nonzero, then note that $\infty$ gets mapped to $\frac{a}{c}$, which needs to be in $\mathbb{R}\cup\{\infty\}$. So either $c=0$ or $\frac{a}{c}$ is real. Well, $a$ is already nonzero real, so $c$ is real too.

Now look at $f^{-1}(z)=\frac{dz-b}{-cz+a}$. We already know $c$ can be assumed real, so the same argument shows $d$ can be assumed real.

Finally your observation about the image of $1$ shows that if $a$, $c$, and $d$ are real, then $b$ is real too.

Edited.

I'm assuming we want $\mathbb{R}\cup\{\infty\}$ to be mapped to $\mathbb{R}\cup\{\infty\}$. The desired conclusion is that we can find $a',b',c',d'\in\mathbb{R}$ such that $$\frac{az+b}{cz+d} = \frac{a'z+b'}{c'z+d'}\quad\text{for all }z\in\mathbb{C}.$$

By plugging in $0$, we conclude that either $d=0$ or $\frac{b}{d}\in\mathbb{R}$.

If $d=0$, plugging in $\infty$ gives $\frac{a}{c}$, hence $\frac{a}{c}=\alpha\in\mathbb{R}$ (or $c=0$, in which case the transformation just gives $z\mapsto \infty$ for all $z$, and we can rewrite it as $\frac{1}{0z+0}$). We can rewrite the transformation as: $$\frac{az+b}{cz} = \alpha + \frac{b}{cz}.$$ Plugging in $z=1$ gives $\alpha+\frac{b}{c}\in\mathbb{R}$, hence $\frac{b}{c}=\beta\in\mathbb{R}$; so we can rewrite $$\frac{az+b}{cz+d} = \frac{az+b}{cz} = \alpha + \frac{\beta}{z} = \frac{\alpha z+\beta}{1z+0},$$ and we are done.

If $b=0$, then composing with $z\mapsto \frac{1}{z}$ we can repeat the argument above. So we may assume that $d\neq 0$ and $b\neq 0$.

Then $\frac{b}{d}=\beta\in\mathbb{R}$, so we can rewrite as $$ \frac{az+b}{cz+d} = \frac{az+\beta d}{cz+d},\quad\beta\in\mathbb{R}.$$ Plugging in $\infty$ we get $\frac{a}{c}=\alpha\in\mathbb{R}$, so we can write $$ \frac{az+b}{cz+d} = \frac{\alpha cz + \beta d}{cz+d}.$$ Plugging in $1$ we get $$\frac{\alpha c + \beta d}{c+d} = \alpha + \frac{(\beta-\alpha)d}{c+d}.$$ Since $\alpha$ and $\beta$ are real, this is a real number (or $\infty$) if and only if $\frac{d}{c+d}$ is a real number (or $\infty$), if and only if $\frac{c+d}{d} = \frac{c}{d}+1$ is a real (or $\infty$), if and only if $\frac{c}{d}$ is a real number (cannot be $\infty$, since $d\neq 0$).

Thus, $c=\gamma d$ with $\gamma\in\mathbb{R}$, so $$ \frac{az+b}{cz+d} = \frac{\alpha cz + \beta d}{cz+d} = \frac{\alpha\gamma dz + \beta d}{\gamma dz + d} = \frac{\alpha\gamma z+ \beta}{\gamma z + 1},$$ with $\alpha,\beta,\gamma\in\mathbb{R}$, as desired.