Most of these identities have a different explanation, via holonomic functions, aka solutions of linear differential equations with polynomial coefficients. And, in the end, it is all a matter of linear algebra and bounding dimensions. This is very well explained in slides by Bruno Salvy (where trigonometrics are explicitly used for the simple examples, before moving on to the 'fun' stuff).
To a certain extent, there is still a lot of geometry still going on in the above, but the geometry of the parameter space, rather than looking directly at the traces of the functions.
Many details of the method can be found in the paper A non-holonomic systems approach to special function identities by Chyzak, Kauers and Salvy.
I asked Mathematica about the boundary behavior. First, the $P^{1/4}_{\nu } $ solution: We have, for $\epsilon \searrow 0$,
$$
(\sin \epsilon )^{1/4} P^{1/4}_{\nu } (\cos \epsilon ) = \frac{2^{1/4} }{\Gamma(3/4)} + O(\epsilon^{2} )
$$
and
\begin{eqnarray*}
(\sin (\pi -\epsilon ))^{1/4} P^{1/4}_{\nu } (\cos (\pi -\epsilon) ) &=& \frac{2^{3/4} \pi }{\Gamma(3/4)\Gamma(-\nu )\Gamma(1+\nu )} \\
& & -\frac{2^{1/4} \pi }{\Gamma(5/4)\Gamma(-1/4-\nu )\Gamma(3/4+\nu )} \sqrt{\epsilon } \\
& & + O(\epsilon^{2} )
\end{eqnarray*}
So, the $P^{1/4}_{\nu } $ solution automatically satisfies the boundary condition at $\theta =0$, whereas at $\theta =\pi $, we have to eliminate the term proportional to $\sqrt{\epsilon } $. That determines the eigenvalues: We need either $-1/4-\nu $ to be a negative integer or 0, or $3/4+\nu $ to be a negative integer or 0.
The specification of $\nu $ in the OP suggests the constraint $\nu \geq -1/2$; this excludes the second alternative, and therefore we obtain the spectrum $\nu = n-1/4$, $n=0,1,2,3,\ldots $.
The $Q^{1/4}_{\nu } $ solution, on the other hand, exhibits behavior proportional to $\sqrt{\epsilon } $ at $\theta=0$, with coefficient
$$
\frac{\pi^{2} }{2^{1/4} } \frac{\cos ((4\nu +1)\pi /8) \Gamma (-\nu /2 -1/8) \Gamma (\nu /2 +9/8) - \sin ((4\nu +1)\pi /8)\Gamma (-\nu /2 +3/8) \Gamma (\nu /2 +5/8)}{\Gamma (5/4) \Gamma (-\nu /2 -1/8) \Gamma (-\nu /2 +3/8) \Gamma (\nu /2 +3/8)\Gamma (\nu /2 +7/8)}
$$
A plot as a function of $\nu $ suggests that, for $\nu \geq -1/2$, this is positive and monotonically rising (I have not attempted to verify this analytically); the behavior proportional to $\sqrt{\epsilon } $ at $\theta=0$ can therefore not be eliminated, nor can it be compensated by admixture of the $P^{1/4}_{\nu } $ solution. Thus, there are no further solutions involving $Q^{1/4}_{\nu } $.
In summary, the complete spectrum is given by $\nu = n-1/4$, $n=0,1,2,3,\ldots $, or, in terms of $\lambda $, $\lambda = n(n+1/2)$, $n=0,1,2,3,\ldots $.
Best Answer
I find it useful to represent the Legendre function in terms of a hypergeometric function, using a formula from Wikipedia,
$$f(\theta)=\frac{ (1+\cos \theta)^{\mu/2} \, _2F_1\left(-\nu,\nu+1;1-\mu;\frac{1}{2} (1-\cos \theta)\right)}{\Gamma (1-\mu)\sin^\mu\theta(1-\cos \theta)^{\mu/2}}.$$
Then Mathematica gives me the small-$\theta$ expansion $$f(\theta)=\frac{ 4 (1-\mu)+\theta^2 \bigl(\tfrac{1}{3}(1-\mu) \mu- \nu (\nu+1)\bigr)+{\cal O}(\theta^4)}{2^{2-\mu}\theta^{2\mu} \Gamma (2-\mu)}.$$ The leading order term is $\nu$-independent, as noted in the OP, the next order term does depend on $\nu$.