When we plot wavefunctions(not its absolute square), we have two options. We can either plot real and imaginary parts separately, or choose a basis so that all basis functions are real-valued (We can do this whenever there is a time reversal symmetry). The latter is what they often do, especially in chemistry books.
For l=1 case, Y_10 is already real, whereas Y_11 and Y_1,-1 are complex.
Y_10 ~ cos(theta) = z/r
Y_11 ~ sin(theta)exp(i*phi) = (x+iy)/r
Y_1,-1 ~ sin(theta)exp(-i*phi) = (x-iy)/r
Yet, if we take the linear combinations of Y_11 and Y_1,-1 we have real valued functions
Y_10 ~ cos(theta) = z/r
Y_11 + Y_1,-1 ~ sin(theta)cos(phi) = x/r
Y_11 - Y_1,-1 ~ sin(theta)sin(phi) = y/r
Now, we have all real orbitals! These orbitals are so called p_z, p_x and p_y orbitals.
The rule of thumb is that you leave m=0 state alone, and take the sum of and difference between all m and -m states for m!=0.
I guess you might want to do the same thing for l=2 and verify that we get x^2-y^2, 3z^2-r^2, xy, yz, zx states.