I see. I can offer a more complicated solution to your problem, but I am definitely interested to hear of a potentially simpler answer if anyone has one.
The MRF you are referring to is the same thing as the terms I mentioned. There’s a very good tutorial on implementing a model with an MRF (aka CAR) component here but it’s not trivial to understand/use if you have no exposure to spatial statistics. The good news is that the MRF model from that link uses some tricks to make sampling much faster than if you tried to construct a Potential term with some ad hoc adjustments.
The challenge is most Bayesian models for gridded data allowing for spatial correlation can’t be specified in terms of sequential correlations (i.e. a directed graphical model) but instead have to be an undirected graphical model requiring definition of the joint likelihood over all sites simultaneously. This is unfortunately more complicated than simply specifying sequential correlations. You can think of the conditional autoregression (CAR) model from the above link as a spatially correlated intercept to be added to the other regression component. The CAR model is also a special case of a Gaussian process (GP) for which PyMC3 has very robust modeling tools.
If you want a quick sketch of what the CAR or GP model might look like then I am happy to show that.