In almost every proof of the mean value formula for harmonic functions (e.g. Evans 2.Theorem 2) one finds the following calculation
$$ \frac{1}{|{\partial B(0,1)}|}\int_{\partial B(0,1)} D u(x+rz)\cdot z\ dS(z) = \frac{1}{|{\partial B(x,r)}|}\int_{\partial B(x,r)} D u(y)\cdot \frac{y-x}{r}\ dS(y)\tag{1},$$
where $u$ is a harmonic function, $r>0$ and $S$ denotes the surface measure.
Obviously, the transformation $x+rz=y$ is used. However, in the literature it is not mentioned why such an application of the transformation formula is allowed. Particularly, the (classical) transformation formula is not applicable to surfaces.
To justify this transformation rigorously, one performs two transformations related to the translation invariance of the Hausdorff/surface measure and uses that $r^{n-1}|\partial B(0,1)|=|\partial B(0,r)|$.
Is there, besides using polar coordinates, another more elementary proof for (1) or am I just thinking too complicated?
You can do this. Let $v_\lambda(x)=u(\lambda x)\,,$ then $v_\lambda$ is also harmonic if $u$ is, so: $\displaystyle \Delta v_\lambda=0\implies\int_{B_r(x)}\Delta v_\lambda \,dV=0\implies \int_{\partial B_r(x)}\frac{\partial v_\lambda}{\partial r} \,r^{n-1}\,dS=0\implies \int_{\partial B_r(x)}\frac{\partial v_\lambda}{\partial \lambda} \,r^{n-1}\,dS=0\implies \int_{\partial B_r(x)}v_\lambda \,dS=\text{constant}$
as a function of $\lambda\,,$ and this last equality implies the mean value theorem (note we used that $r$ is constant on spheres to cancel it out). Do you have the same problem with it?