We analyze several types of Galerkin approximations of a Gaussian random field mathscr {Z}:mathscr {D}times varOmega rightarrow mathbb {R} indexed by a Euclidean domain mathscr {D}subset mathbb {R}^d whose covariance structure is determined by a negative fractional power L^{-2beta } of a second-order elliptic differential operator L:= -nabla cdot (Anabla ) + kappa ^2. Under minimal assumptions on the domain mathscr {D}, the coefficients A:mathscr {D}rightarrow mathbb {R}^{dtimes d}, kappa :mathscr {D}rightarrow mathbb {R}, and the fractional exponent beta >0, we prove convergence in L_q(varOmega ; H^sigma (mathscr {D})) and in L_q(varOmega ; C^delta (overline{mathscr {D}})) at (essentially) optimal rates for (1) spectral Galerkin methods and (2) finite element approximations. Specifically, our analysis is solely based on H^{1+alpha }(mathscr {D})-regularity of the differential operator L, where 0<alpha le 1. For this setting, we furthermore provide rigorous estimates for the error in the covariance function of these approximations in L_{infty }(mathscr {D}times mathscr {D}) and in the mixed Sobolev space H^{sigma ,sigma }(mathscr {D}times mathscr {D}), showing convergence which is more than twice as fast compared to the corresponding L_q(varOmega ; H^sigma (mathscr {D}))-rate. We perform several numerical experiments which validate our theoretical results for (a) the original Whittle–Matérn class, where Aequiv mathrm {Id}_{mathbb {R}^d} and kappa equiv {text {const.}}, and (b) an example of anisotropic, non-stationary Gaussian random fields in d=2 dimensions, where A:mathscr {D}rightarrow mathbb {R}^{2times 2} and kappa :mathscr {D}rightarrow mathbb {R} are spatially varying.