Wondering if there is a high-school level approach that doesn't involve manipulation of special functions or convolutions or complex analysis.
It is easy to prove this holds for positive integer values of $r$, but trying to extend this beyond seems difficult.
I thought of maybe using polynomial properties (vanishing at infinitely many points implies identically 0) but I am not sure if this is even a workable proof.
We can rewrite the product as
$$\int_0^{\frac{\pi}{2}}\sin^r\theta\:d\theta\int_0^{\frac{\pi}{2}}\sin^{r-1}\varphi\:d\varphi = \int_0^{\frac{\pi}{2}}\int_0^{\frac{\pi}{2}}(\sin\theta\sin\varphi)^{r-1}\sin\theta \:d\theta \:d\varphi = \iint\limits_{S}y^{r-1}\:dS$$
where $S$ is the portion of the unit sphere in the first octant. By symmetry, the surface integral is the same even if we exchanged $y$ for any of the other coordinate variables thus we have
$$\iint_S y^{r-1}\:dS = \iint_Sz^{r-1}\:dS = \int_0^{\frac{\pi}{2}}\int_0^{\frac{\pi}{2}}\cos^{r-1}\theta\sin\theta \:d\theta\:d\varphi $$
$$=\int_0^{\frac{\pi}{2}}\:d\varphi\int_0^1 t^{r-1}\:dt = \frac{\pi}{2r}$$