Hello—
Just to add another option, there is KernelMatrices.jl, which provides very scalable methods for maximum likelihood estimation. It’s pretty flexible in the sense that you really just need to specify the covariance matrix to the point where individual elements can be computed (and the derivatives if you want a scalable gradient, and the second derivatives if you want a scalable Hessian). There are two example files that demonstrate the pretty small amount of boilerplate code required to use the package for estimation.
I only implemented it for mean-zero fields, but at least in spatial statistics it is not uncommon to try and remove a mean term and then separately estimate the covariance structure. If people really want parametric mean functions included as well I don’t think that it would be too much work.