Uniqueness of the Transformation
Some understanding of group actions can help us show that such a transformation (using real coefficients) must be unique (assume $w\ne z$). If $f$ and $g$ are two such transformations, then $f^{-1}\circ g$ is a transformation that fixes both $z$ and $w$. But writing out the equation for $z$ to be a fixed point of $\frac{az+b}{cz+d}$ reveals (again, for real values $a,b,c,d$) there can only be one fixed point $z$ in the upper half-plane, unless it is the identity. But if it is the identity, $f^{-1}\circ g=\mathrm{id}$ implies $f=g$ are the same transformation.
Geometry of the Transformation
The upper half-plane is a model for hyperbolic geometry, which has a metric. The geodesics are semicircles with diameters on the real axis. Around every point there is a group of "hyperbolic rotations" which fix it and whose derivative there rotates tangent vectors. In particular, around $i$, it is $\mathrm{SO}(2)$ (however the kernel of its action is $\{\pm I_2\}$ so maybe we want to think of it as $\mathrm{PSO}(2)$). Given two points $z$ and $w$ we may construct the unique semicircle arc between them. On that arc is a "hyperbolic midpoint" which is the same hyperbolic distance from both $z$ and $w$ (I expect this is not the midpoint according to Euclidean arclength, though). We can swap $z$ and $w$ by applying a $180^{\circ}$ hyperbolic rotation around this midpoint. (Note that if the midpoint was e.g. $i$, we'd need a $90^{\circ}$ rotation matrix in $\mathrm{SO}(2)$ to achieve this effect.)
Formula for the Transformation
The equations that say $f$ swaps $w$ and $z$ are given by
$$ \frac{aw+b}{cw+d}=z, \qquad \frac{az+b}{cz+d}=w. $$
Clearing denominators yields
$$ \begin{cases}
aw+b = c(wz)+dz \\
az\,+\,b = c(zw)+dw
\end{cases} $$
Subtracting gives $a(w-z)=d(z-w)$ which implies $d=-a$, and the equation becomes
$$ a(w+z)+b=c(wz). \tag{$\ast$} $$
Taking imaginary parts yields $a\mathrm{Im}(w+z)=c\,\mathrm{Im}(wz)$. This allows us to solve:
$$ \begin{cases}
a = \lambda \mathrm{Im}(wz) \\
c = \lambda \mathrm{Im}(w+z)
\end{cases} $$
where $\lambda\in\mathbb{R}^{\times}$. We may solve for $b$ by taking real parts of $(\ast)$ to finally get
$$ \begin{pmatrix} a & b \\ c & d \end{pmatrix}
= \lambda
\begin{pmatrix} \mathrm{Im}(wz) & \mathrm{Re}(wz)\mathrm{Im}(w+z)-\mathrm{Im}(wz)\mathrm{Re}(w+z) \\
\mathrm{Im}(w+z) & -\mathrm{Im}(wz) \end{pmatrix} $$
Special Example
Taking $z=w=i$ gives a $90^{\circ}$ rotation matrix, yielding a $180^{\circ}$ hyperbolic rotation around $i$. This is the limiting case of $z,w\to i$ with midpoint $\to i$.