  Large or very large spatial (and spatio-temporal) datasets have become common
place in many environmental and climate studies. These data are often collected
in non-Euclidean spaces (such as the planet Earth) and they often present
non-stationary anisotropies. This paper proposes a generic approach to model
Gaussian Random Fields (GRFs) on compact Riemannian manifolds that bridges the
gap between existing works on non-stationary GRFs and random fields on
manifolds. This approach can be applied to any smooth compact manifolds, and in
particular to any compact surface. By defining a Riemannian metric that
accounts for the preferential directions of correlation, our approach yields an
interpretation of the ''local anisotropies'' as resulting from ''local''
deformations of the domain. We provide scalable algorithms for the estimation
of the parameters and for optimal prediction by kriging and simulation able to
tackle very large grids. Stationary and non-stationary illustrations are
provided.
