Is there an underlying explanation for the magical powers of the Schwarzian derivative? Given a function $f(z)$ on the complex plane, define the Schwarzian derivative $S(f)$ to be the function
$S(f) = \frac{f'''}{f'} - \frac{3}{2} (\frac{f''}{f'})^2$
Here is a somewhat more conceptual definition, which justifies the terminology.  Define $[f, z, \epsilon]$ to be the cross ratio $[f(z), f(z + \epsilon); f(z + 2\epsilon), f(z + 3\epsilon)]$, and let $[z, \epsilon]$ denote the cross ratio $[z, z + \epsilon, z + 2\epsilon, z + 3\epsilon]$ (in fact this is just 4, but this notation makes my point clearer).  One can ask if $[f, z, \epsilon]$ is well approximated by $[z, \epsilon]$; indeed, it turns out that the error is $o(\epsilon)$.  So one pursues the second order error term and finds that $[f, z, \epsilon] = [z, \epsilon] - 2 S(f)(z) \epsilon^2 + o(\epsilon^2)$.  So the Schwarzian derivative measures the infinitesimal change in cross ratio caused by $f$.  In particular, $S(f)$ is identically zero precisely for Möbius transformations. 
That's all background.  From what I have said so far, the Schwarzian derivative is at best a curiosity.  What is not obvious at first glance is the fact that the Schwarzian derivative has magical powers.  Here are some examples:
First magical power: The Schwarzian derivative is deeply relevant to one dimensional dynamics, stemming from the fact that it behaves in a specific way under compositions.  For example, if $f$ is a smooth function from the unit interval to itself with negative Schwarzian derivative and $n$ critical points, then it has at most n+2 attracting periodic orbits.
Second magical power: It says something profound about the solutions to the Sturm-Liouville equation, $f''(z) + u(z) f(z) = 0$.  If $f_1$ and $f_2$ are two linearly independent solutions, then the ratio $g(z) = f_1(z)/f_2(z)$ satisfies $S(g) = 2u$.
Third magical power: The Schwarzian derivative is the unique projectively invariant 1-cocyle for the diffeomorphism group of $\mathbb{R}P^1$.  This is probably just a restatement of the conceptual definition I gave above, but I'm not sure; in any event, this gives the Schwarzian derivative a great deal of relevance to conformal field theory (or so I'm told).  
I'm sure there are more.  I'm wondering if all of these powers can be explained by some underlying geometric principle.  They all seem vaguely relevant to each other, but the first power in particular seems very hard to relate to the definition in any obvious way.  Does anybody have any insights?
 A: The velocity for a solution $f(x-c\cdot t)$ of the one dimensional heat equation
$$\frac{\partial}{\partial t}f(x-ct)=-c\frac{\partial}{\partial x}f(x-ct)=\alpha\frac{\partial^{2}}{\partial x^{2}}f(x-ct)\;$$
is determined by $$c=\frac{-\alpha D_{x}^{2}\,f(x)}{D_{x\,}f(x)}=-\alpha D_{x}\ln[D_{x}f(x))]\,.$$
Similarly, the velocity for a soliton solution $u(x-ct)$
of the KdV equation
$$u_{t}-u\cdot u_{x}+\frac{1}{12e_{2}}\,u_{xxx}=0\,$$
is determined by
$$c=\frac{D_{x}^{3}u(x)-6e_{2}D_{x}(u(x))^{2}}{12e_{2}D_{x}u(x)}=\frac{1}{6e_2}S_{x}\{h^{-1}(x)\}(g(x))^{2}=-\frac{1}{6e_2}S_{x}\{h(x)\}$$
$$=\frac{e_{2}}{3}\left(\frac{\omega_{2}-\omega_{1}}{2}\right)^{2},$$
where $\,u(x)=h'(x)$ and $\,g(x)=e_{2}(x-\omega_{1})(x-\omega_{2})=\frac{1}{d[h^{-1}(x)]/dx}=h^{'}(h^{-1}(x))$.
Note:
$$g(h(x)) = e_{2} \; (h(x)-\omega_{1})\;(h(x)-\omega_{2}) = h'(x)$$
is a Riccati equation, and the Schwarzian can be expressed as
$$S_{x}\{h(x)\}=D_{x}^{2}\ln[D_{x}h(x)]-\frac{1}{2}\left[D_x\ln[D_{x}h(x)]\right]^{2}=D_{x}^{2}\ln[u(x)]-\frac{1}{2}\left[D_x\ln(u(x))\right]^{2} \; .$$
(See Old and New on the Schwarzian Derivative by Osgood for an excellent survey on the derivative and Peter Michor's papers on the relationships among solutions of the KdV eqn., geodesics of the Virasoro-Bott group, the Schwarzian, and Moebius transformations--the homogeneous solns.)
(Explicit note on the Riccati equation added July 2, 2021,)
A: Let G(k,n) be the k'th order jet group in n variables which consists of the set of k-jets of local diffeomorphisms of R(n) fixing the origin under the operation of composition. In coordinates, this group operation can be written explicitly using the chain rule.
Now consider G(3,1) and G(2,1) and the obvious projection homomorphism from G(3,1) onto G(2,1) induced by jets. This projection splits, that is, there is an injective homomorphism which imbeds G(2,1) into G(3,1). Identifying the image of this injection with G(2,1), we can form the left (say) coset space G(3,1)/G(2,1). Now the expression for the Schwarzian derivative defines "coordinates" on G(3,1)/G(2,1). The details of these computations can be found on pages 152-153 of the book "An Alternative Approach to Lie Groups and Geometric Structures". This construction generalizes to arbitrary dimensions and is a very special case of defining geometric structures as explained on pages 174-182 of this book. 
To summarize here, these splittings arise from homogeneous spaces and define "geometric connections" which are not necessarily connections in the classical sense (keeping in mind that connections can be defined on general principal and vector bundles and are essentially topological objects). Furthermore, these splittings are built into the definition of a geometric structure and therefore there is no need to search for a "special connection" suitable for some geometric structure. Consequently, the magical powers of the Scwarzian derivative is a special case of the magical powers of "connections". 
As an interesting detail, the obvious left action of G(3,1) on G(3,1)/G(2,1) gives the transformation rule of the Schwarzian derivative found in classical textbooks. In affine case (the action of G(2,n) on G(2,n)/G(1,n)), for instance, we get the well known transformation rule of the "connection components" on the tangent bundle which is historically the starting point of the theory of connections on vector bundles! The action of G(1,n) on G(1,n)/G(0,n) = G(1,n) gives absolute parallelism studied in detail in this book.
A: One very simple interpretation is possible, if I am not mistaken: 
Characterizing conformal nets, or more general harmonic nets having in mind
the lift to minimal surfaces, it turns out that
there exists a simple uniform relationship of four fundamental entities up
to normalization for orthogonal trajectories:


*

*The sum of the changes of geodesic curvature is zero.

*The difference of the changes of geodesic curvature is the real
part of the Schwarzian derivative.

*The difference of the changes of geodesic acceleration is the
imaginary part of the Schwarzian derivative.

*The sum of the changes of geodesic acceleration is the Gaussian
curvature.
http://vixra.org/abs/1301.0071
A: The following book may contain insights of the type you are looking for, at least for the second and third interpretation:

V. Ovsienko, S. Tabachnikov, Projective differential geometry old and new.
  From the Schwarzian derivative to the cohomology of diffeomorphism groups. Cambridge Tracts in Mathematics, 165. Cambridge University Press, Cambridge, 2005. ISBN: 0-521-83186-5 MR

Although I didn't read the book, earlier articles by the same authors, including a popular exposition in the now defunct "Quant" magazine, were amazingly insightful. Here is how the authors begin their story:

Every working mathematician has encountered the Schwarzian derivative at
  some point of his education and, most likely, tried to forget this rather scary
  expression right away. One of the goals of this book is to convince the reader
  that the Schwarzian derivative is neither complicated nor exotic, in fact, this
  is a beautiful and natural geometrical object.


I'd like to complement the "third magical power" with the statement that the Schwarzian derivative is a cocycle of the diffeomorphism group of the circle $S^1$ with values in the quadratic differentials on $S^1$, and since the latter space can be identified with the dual space of the vector fields on $S^1,$ it encodes a central extension of $\text{Diff}(S^1);$ the infinitesemial version of this central extension is the Virasoro algebra.
A: Nobody pointed out What is ...  Schwarzian Derivative? (Notices of AMS Jan 2009), which succinctly explains quite a lot.
A: It's a cocycle in group cohomology.
$S[f\circ g]= S[f]^g + S[g]$
which is the cocycle condition! The group is Aut(meromorphic functions) and the other group is meromorphic functions.
A: Like may people (but not all people), I have trouble thinking in terms of formulas such as that for the Schwarzian. For me, a geometric image works much better. I'll describe a geometric picture, similar to what I discussed in my paper "Zippers and univalent functions" (and can be found elsewhere as well, but I don't have a good sense of references).
The group of Moebius transformations is 3 dimensional, and for any locally-defined diffeomorphism f of R or any locally-defined holomorphic map of C, you can fit the value, the first derivative and the 2nd derivative at any point by a unique Moebius transformation, the osculating Moebius transformation at the point.   From this, you can make a recipe to extend the diffeomorphism into the upper half plane or upper half space models for hyperbolic geometry: map each vertical line according to the osculating Moebius transformation at its base.  When you do this, vertical lines are mapped isometrically, but the metric is necessarily distorted in the horizontal directions unless f is a Moebius transformation. The Schwarzian derivative gives the asymptotic behavior of this distortion.  For real maps, if the Schwarzian is negative, the hyperbolic rays are bent away from each other.  A hyperbolic line perpendicular to the vertical lines is mapped to a curve that (in terms of the hyperbolic metric) bends downward.  If you consider any interval on R, it has a natural projectively-invariant metric, called the Hilbert metric, and identified with the 1-dimensional hyperbolic metric in this case. The bending  implies that the metric is expanded by f, relative to the Hilbert metric of its image.  For example, $\log(t)$ is an arc-length parametrization for $(0,\infty)$.  The map $x \rightarrow x^k$ expands the parameter by a factor of $k$.   The expanding property is highly significant for analyzing dynamics.
In the complex analytic case, the Schwarzian is a holomorphic quadratic differential, and can be geometrically indicated by two perpendicular foliations:  one set of streamlines where the quadratic form takes positive real values, and a perpendicular set of lines where it takes negative real values. Whenever the quadratic form has a simple zero, the foliations have singularities where the lines make a Y pattern, branching in 3 ways.  At any critical point of f, there is a double pole of the Schwarzian, where the positive real foliation circles around and the negative real foliation is radial.  These lines show how the extension of f bends surfaces asymptotically near the complex plane; if you start with an umbilic surface such as a plane, horosphere or equidistant surface, they show the asymptotic pattern for lines of curvature for the image of the surfaces via the extension of f.  For example, you can visualize z -> z^2 as extending by mapping the hyperbolic cylinders around the z axis (they appear as cones in upper half space), wrapping them twice around the vertical axis (thus stretching the circumference by a factor of 2) , and stretching vertically by a factor of 2.  The curvature along the meridians is increased, and the curvature in directions parallel to the axis is decreased.  Of course, the Poincare metric of a small disk is preserved, but you can still see the metric effect in terms of the nonconformality of the behavior on surfaces in hyperbolic space.  You can also see it by the shape of the image of a small circle on C.  To 2nd order, it remains round, but there is a 3rd order effect that makes it elliptical, where the short axis is the direction in which the Schwarzian is negative.   
When you look at computer plots of the quadratic differentials for holomorphic maps, they pop into 3 dimensions, strongly suggesting the geometry of some families of surfaces that can be associated to a holomorphic map of C.
 (source: Wayback Machine)
The Schwarzian for the rational function $f(z)(= (z^3-3 z -1)/(z^3+1)$. It's challenging to make a revealing plot for a complex rational function like $f$, since it maps 3 times over the Riemann sphere, but the Schwarzian is easy to draw, and shows how an extension of f bends hyperbolic space. The critical points of f(z) are surrounded by circular positive-real circles where the Schwarzian has a double pole, indicating how they are wrapped twice around a core singularity. A typical zero is visible in the center, where the bending branches three ways. The zeros and poles of the function f itself are not visible, since $0$ and $\infty$ have no special significance in the geometry of $S^2 = CP^1$.
Enough ... there are endless mysteries to the Schwarzian, but these geometric images are helpful to me.   I'll second Victor's suggestion: I haven't read the Osvienko- Tabachnikov book either, but I'm confident from prior experience that it's interesting material.
Addendum. Since you expressed interest in the real case, here's one way to illustrate it. The idea is simple: take a standard family of horocycles in the domain, in this case circle of constant height tangent to the real line, and push them forward by the osculating Moebius transformation. The image is actually determined by just $f$ and $f'$, but to calculate (rather than just see) the envelope would require $f''$. This picture is for $x^3 - 3 x$ (which folds $[-2,2]$ 3 times over its image), in the interval $[-2.1,2.1]$  The fat shape with downard hyperbolic curvature of the envelope in contrast to upward curvature of the envelope in the domain, suggests a caterpillar outgrowing its skin and demonstrates the negative Schwarzian.
 (source: Wayback Machine)
A: The Schwarzian derivative encodes the adjoint action of the Bott-Virasoro group. One version of it is 
${Diff}\_{\mathcal S}(\mathbb R)\times \mathbb R$
(here $\mathcal S$ stands for "rapidly falling towards the identity") with multiplication
$$
\binom{\phi}{\alpha}.\binom{\psi}{\beta} = \binom{\phi\circ\psi}{\alpha+\beta+c(\phi,\psi)}
$$
where the Bott cocycle is: 
$$
c(\phi,\psi) = \frac12\int\_{\mathbb R} \log(\phi'\circ \psi)\,d\log(\psi').
$$
This is alluded to on page 22ff of the book of Ovsienko and Tabashnikov mentioned in Viktor Protsaks answer. A short exposition along the lines of this answer is on page 55 of 
(here). I think that this encodes much of the magic of the Schwartzian derivative. Dualize it to the coadjoint action and note that many coadjoint orbits are of the form $Diff\_{\mathcal S}(\mathbb R)/PSL(2,\mathbb R)$ to see the projective properties. Etc.
A: Let $f:\mathbf{R}\mapsto\mathbf{P}^1$ be a function realised as $y=f(x)$, where $y$ is an affine coordinate on $\mathbf{P}^1$. Since the target space is the projective line it is desirable to lift $f$ to homogeneous coordinates, i.e. to choose $\mathbf{v}=(u(x), v(x))$ such that $u/v=f$. (Maybe it is not quite right to call this homogenous coordinates, rather this is a lifting of $f$ to the line bundle over $\mathbf{P}^1$, i.e. converting $f$ to a parameterised curve in $\mathbf{R}^2$.) Differentiating the constraint gives $$\frac{u'v-uv'}{v^2} = f',$$ which suggests taking $u'v-uv'=1$ as an additional constraint (assuming $f'>0$). This fixes $u$ and $v$ to be $$u=\frac{f}{\sqrt{f'}}$$ $$v=\frac{1}{\sqrt{f'}}$$
The added constraint says that the angular momentum of the lifted curve is constant, $\mathbf{v}'\times\mathbf{v}=1$ and this means that the acceleration vector is parallel to the position vector, $\mathbf{v''}=\lambda\mathbf{v}$. The proportionality factor $\lambda = \lambda(x)$ is $$\lambda = \frac{u''}{u}=\frac{v''}{v}=-\frac{1}{2}\frac{f'''f'-\frac{3}{2}(f'')^2}{(f')^2},$$ which is the Schwarzian with an extra factor of $-1/2$.
The proportionality factor $\lambda$ is invariant under linear transformations of $\mathbf{R}^2$. Since a projective coordinate change on the target space lifts to such a linear transformation, the Schwarzian is independent of the choice of the affine coordinate $y$ and is a projective invariant.
