Let $\Bbb S^{n-1}$ denotes the unit sphere in $\Bbb R^n$ and let $O(n)$ be the group of all orthogonal linear transformations on $\Bbb R^n$. For any (outer) measure $\mu$ on $\Bbb R^n$, the pushforward of $\mu$ by $T\in O(n)$ is defined as $T_{\#}\mu(E) := \mu(T^{-1}(E))$. It is clear that the pushforward under $T\in O(n)$ maps a probability measure on $\Bbb S^{n-1}$ to another probability measure on $\Bbb S^{n-1}$.
Suppose that $f\in C(\Bbb S^{n-1})$ satisfies the condition $$ \int_{\Bbb S^{n-1}} f \,d\mu = \int_{\Bbb S^{n-1}} f \,dT_{\#}\mu \quad \text{for all } \ T\in O(n), \tag{*}\label{*} $$ is it true that $$ \int_{\Bbb S^{n-1}} f \,d\mu = \frac{1}{\mathcal H^{n-1}(\Bbb S^{n-1})} \int_{\Bbb S^{n-1}} f \,d\mathcal H^{n-1} =: \bar f, \tag{**}\label{**} $$ where $\mathcal H^{n-1}$ is the $n-1$ dimensional Hausdorff measure?
Of course, the case where $\mu$ is a Dirac mass $\delta_x$ for some $x\in\Bbb S^{n-1}$ is clear, and, in fact, we can trivially conclude that $f$ must be a constant function. One might be tempted to conjecture that the condition $\eqref{*}$ implies that $f$ must be a constant function, but that is false. We can already find a counterexample to the stronger conjecture by considering $\mu$ to be a sum of Dirac masses: $\mu = \frac1n (\delta_{e_1} + \cdots + \delta_{e_n})$ and $f(x) = x\cdot Ax$ for any $n\times n$ matrix $A$. Indeed, in this case we have $$ \int_{\Bbb S^{n-1}} f \,d\mu = \frac 1n \sum_{i=1}^n e_i\cdot Ae_i = \frac{\text{Tr}(A)}{n}, $$ and the statement that $\eqref{*}$ holds follows from the well-known fact that the trace of a matrix $A$ is invariant under a change of an orthogonal basis. Clearly, $f$ is not a constant function, but it is not hard to show that $$ \frac{\text{Tr}(A)}{n} = \frac{1}{\mathcal H^{n-1}(\Bbb S^{n-1})} \int_{\Bbb S^{n-1}} x\cdot Ax \,d\mathcal H^{n-1}(x) $$ using the divergence theorem and that fact that $\text{Tr}(A) = \text{div}\, Ax$.
This question is largely inspired by the example regarding the trace that I mentioned above. I restrict $f$ in this question to be a continuous function since I feel like there might be a beautiful proof using the duality between continuous functions and Radon measures on $\Bbb S^{n-1}$.
I am aware of the existence of the things called measures on homogeneous spaces that generalize Haar measures, but aside from that, I know next to nothing about it. I welcome a proof that uses ideas from that theory, but would like to read a proof that doesn't use it as well, even if it might be longer and less elegant.
Edit: Take a look at BigbearZzz' answer for a more detailed explanation!
Since the integration of $\int_{S^{n-1}} f d \mu$ is invariant under pushforwards by $T$, we can average over all $T \in O(n)$ via the (normalized) Haar measure $\nu$ on $O(n)$ and interchange the integrals by Fubini's theorem. \begin{align}\int_{S^{n-1}} f d \mu &= \int_{O(n)} \int_{S^{n-1}} f(x) \, d(T_{\#} \mu)(x) \, d \nu(T) = \int_{O(n)} \int_{S^{n-1}} f(T(x)) \, d \mu(x) \, d \nu(T) \\ &= \int_{S^{n-1}} \int_{O(n)}f(T(x)) \, d \nu(T) \, d \mu(x) = \int_{S^{n-1}} \int_{O(n)}f(e_x(T)) \, d \nu(T) \, d \mu(x) \\ &= \int_{S^{n-1}} \def\avint{\mathop{\,\rlap{-}\!\!\int}\nolimits} \avint_{S^{n-1}} f(y) \, d \mathcal{H}^{n-1}(y) \, d \mu(x) = \avint_{S^{n-1}} f(y) \, d \mathcal{H}^{n-1}(y) \\ &= \bar{f}. \end{align} Here, the fourth equality holds since $x$ is fixed, so integrating over the Haar measure on $O(n)$ really averages over the sphere.
More precisely, $O(n)$ is the isometry group of $S^{n-1}$ (for example by this post), and so for arbitrary fixed $x \in S^{n-1}$ and evaluation map $$e_x : O(n) \to S^{n-1}, \quad T \mapsto e_x(T) := T(x),$$ we have $$ (e_x)_{\#} \nu = \frac{\mathcal{H}^{n-1}}{\mathcal{H}^{n-1}(S^{n-1})}.$$ (For example by this post).