Chemists use the linear combinations of hydrogen atom wavefunctions which are real because this is simply much more convenient to think about. Notice that for hydrogen, all the orbitals for a given $n$ are degenerate, so taking these linear combinations does not affect the energy, and thus thinking about the shapes of these orbitals is potentially more useful than thinking about the shape of something which is complex.
You are correct that taking the real-valued wavefunctions causes them to no longer be eigenfunctions of $L_z$. This is a bit frustrating, but in practice we never use these orbitals for anything, so it does not really matter. I used to think organic chemists took the shape of orbitals a bit too seriously (and they maybe do), but the notion of the shape of an orbital is actually recovered in a more rigorous manner.
That is, molecular orbital theory is just a qualitative form of Hartree-Fock (HF) theory. In HF, the Fock matrix is unitarily invariant, so we are free to transform our orbitals by any unitary rotation. These orbitals can be written as functions in $\mathbb{R}^3$ as desired. Thus, there exists schemes by which molecules orbitals can be "localized" and they often take the shape of orbitals which look a whole lot like the real-valued hydrogen atom orbitals you mention.
Thus, the use of these real-valued orbitals as a notion of shape is justified after the fact from HF theory and beyond.
In practice, this never really matters if one is being rigorous because any attempt at saying exactly which orbitals are involved in a bond is thwarted by unitarily transforming those orbitals. Thus, many chemists use the shape of an orbital as a heuristic because thus far, it has proven to be a very useful conceptual tool. I think pretty much all chemists know this is not completely correct, but neither is anything in physics.
Basically, the shape of orbitals provide a useful model.
