Kronecker Structured Covariances may help for GP with spatial data. To use LatentKron as prior, have a look at this pymc-example.
Also, a recent blog post “Modeling spatial data with Gaussian processes in PyMC” from from PyMC Labs Blog may be useful as well.