I'll try to expand on my comment, and address some apparent misconceptions in your question.
Lets start with your comment 'we should have been better off plotting $E-r$' -
why should that be? When we look at electron orbits around an isolate nucleus,
we get discrete wavefunctions (the s, p, d, f, ... functions). While the probability
distribution can be looked at as a function of the electron-nucleus distance $r$, the
energy of the orbital is a function solely of the quantum numbers. It is not a function
of $r$ at all. We also look for wavefunctions that are centered on the nucleus - there
is no consideration for translation of the nucleus in atomic energy levels, so there is no linear momentum in an atomic orbital.
Now lets go to a crystalline solid. True, we abandon atomic orbital
descriptions of electron energy levels. Why? Not just 'because', we do so since we
have a description of a solid as a Bravais lattice, a space-filling description using
lattice vectors. We look for electron wavefunctions that are Bloch functions which
take advantage of this representation and insist that the Bloch functions be constrained
by the lattice. That is, the Bloch functions are invarient under translation by arbitrary
linear combinations of the unit vectors. The result of solving this (non-atomic) wave equation
is a wavefunction that explicitly has linear momentum in it. A given energy level has
a linear momentum associated with it since translation with the lattice is allowed (required).
Further, since the Bloch function is periodic in the
lattice spacing $a$ (a requirement to be a Bloch function), the result is also restrictions on
momentum, or, more specifically, that a Bloch function with some momentum $k$ is the same as if it had
a momentum $k + {{2\pi}\over{a}}$ - the translational symmetry implies a symmetry in momentum space.
Digressing a bit, we are quite comfortable dealing with $E$ vs $k$ classically,
in that $E = {1\over 2} mv^{2}$. Since classical momentum is $p=mv$, we get
$E = {p^{2} \over 2m}$. Thus, a 'free' classical electron has an $E$ vs $k$ diagram that
is a parabola.
Now, lets put a 'free' electron in a crystal. We get a nice parabola in $E$ vs $k$, since being
a free electron means it is classical in nature. But, wait, the classical 'free' electron
is not a Bloch function - it has no relationship with the lattice. So, lets make in 'nearly'
free, forcing it to be a Bloch function with little other impact. Again, you get a parabola in $E$ vs $k$.
OK, good. Now apply the momentum symmetry from above - take that parabola centered on $k = 0$, and
translate it by multiples of $\pm {{2\pi}\over{a}}$. Those are all equivalent levels. The result is you
get a bunch of parabolas intersecting at points at the zone edges.
Next, lets look at what 'nearly' free really means. It means there is some small interactions
going on. So, focus in on the intersections in $E$ vs $k$ space - you have two solutions
trying to be in the same place for an electron, which is a fermion. This is a problem. You now have some
interaction (the 'nearly' part) which will lead to level splitting - you can't actually have
the parabolas cross. You get a forbidden zone, the width of which depends on the strength of the
interaction.
Going to even-less-free electrons (which is most crystals) leads the the wide range of band structures
that you can find in nature, including having minima of bands away from the symmetry points (like in
silicon), minima of an upper band below the maxima of a lower band but separated in $k$ (semimetals like tin),
and so forth and so on.