Gaussian random fields play an important role in many areas of science and engineering. In practice, they are often simulated by sampling from a high-dimensional multivariate normal distribution, which arises from the discretisation of a suitable precision operator. Existing methods such as Cholesky factorization and Gibbs sampling become prohibitively expensive on fine meshes due to their high computational cost. In this work, we revisit the Multigrid Monte Carlo (MGMC) algorithm developed by Goodman & Sokal (Physical Review D 40.6, 1989) in the quantum physics context. While the authors of this paper conclude that MGMC does not overcome critical slowing down in simulations of field theories near phase transitions, we demonstrate here that it has the potential to significantly accelerate sampling in spatial statistics. The class of Gaussian Random Fields we consider includes those with Matérn covariance, but is more general in that it also allows for non-stationary covariance functions. To show that MGMC can overcome the limitation of existing methods, we establish a grid-size-independent convergence theory based on the link between linear solvers and samplers for multivariate normal distributions, drawing on standard multigrid convergence arguments. We then apply this theory to linear Bayesian inverse problems. This application is achieved by extending the standard multigrid theory to operators with a low-rank perturbation. Moreover, we develop a novel bespoke random smoother which takes care of the low-rank updates that arise in constructing posterior moments. In particular, we prove that Multigrid Monte Carlo is algorithmically optimal in the limit of the grid-size going to zero. Numerical results support our theory, demonstrating that Multigrid Monte Carlo can be significantly more efficient than alternative methods when applied in a Bayesian setting.