Usually the procedure is the other way around, i.e. one starts with the observables and then (if there are sufficiently many commuting observable) one constructs eigenstates of the operators corresponding to the observables. A wavefunction will then be some combination of the eigenstates specified by the initial conditions of the problem.
Since the eigenvalues of the operators play the role of labels for the eigenstates, the choice of observables is often dictated by the physics of the problem. The energy is an important quantity; the corresponding operator is the Hamiltonian so one usually constructs states that are eigenstates of the Hamiltonian. For 1d problem this is enough.
For 2d or 3d problems one must look for additional observables. In the case of a 3d spherically symmetric potential, those additional observables are the square of the angular momentum operator $\vec L\cdot \vec L$ (with operator $\hat L^2$) and one component of the angular momentum, usually $L_z$ (with operator $\hat L_z$). Any other component of $\vec L$ would work but spherical coordinates "single out" the $\hat z$ component as having a particularly simple expression $L_z=L\cos\theta$ so this choice is the more common. In more mathematical terms:
\begin{align}
\hat H\psi_{n\ell m}(r,\theta,\phi)&=E_{n,\ell}\psi_{n\ell m}(r,\theta,\phi)\, , \\
\hat L^2\psi_{n\ell m}(r,\theta,\phi)&=\hbar^2\ell(\ell+1)\psi_{n\ell m}(r,\theta,\phi)\, , \\
\hat L_z\psi_{n\ell m}(r,\theta,\phi)&=\hbar m \psi_{n\ell m}(r,\theta,\phi)\, .
\end{align}
Note that $E, \hbar^2\ell(\ell+1)$ and $\hbar m$ are number (eigenvalues in fact), not observables.
There are of course other observables in your system: $L_y$ and $L_x$ are obvious, but one can think of various other combinations of $p_k$'s and $q_j$'s. However, in general, for observables ${\cal O}$ other than $E$, $\vec L\cdot\vec L$ or $L_z$
$$
\hat {\cal O}\psi_{n\ell m}(r,\theta,\phi)\ne o\,\psi_{n\ell m}(r,\theta,\phi)\, , \tag{1}
$$
but it still makes sense to compute average values and matrix elements, such as
$$
\langle {\cal O}\rangle = \int dV \psi^*_{n\ell m}(r,\theta,\phi)\hat{\cal O}\psi_{n\ell m}(r,\theta,\phi)
$$
(There are exceptions to (1) if you can find observables that commute with the above three, but this does not always happen.)
Nevertheless, you can verify that in systems with spherically symmetric potentials, the basis functions are usually eigenstates of the energy, the magnitude of the angular momentum, and the $z$-component of angular momentum. These three observables commute and are enough to completely specify an eigenstate. Moreover, they form a complete set of commuting observables, meaning that any wavefunction can be expanded as a linear combinations of eigenstates of these operators.
If you choose to include spin, then this set of observable is no longer complete, and you need to include spin operators to completely specify eigenstates.
