Everything you’re saying sounds doable, you just need to replace “loop over” with “compute over a vector” in your planning. As long as there’s no inter-temporal dependencies that require recursive computation, you can treat the time dimension as just another index value and do the vector computation over the whole dataframe in one shot. If you’re not familiar with vectorized computation with numpy, it would be worthwhile to check out some resources on that, because pytensor (the pymc computational backend) works exactly the same way. If you can write down your model in numpy (without loops) it’s a basically a one-to-one mapping to a PyMC model, and you can drop in random variables wherever you need.
So C(x,y,z) = Q \cdot cc + \epsilon? In that case you can use C(x,y,z) to parameterize some distribution, for example as the mean of the normal (if \epsilon is mean zero normal)
Do you mean that you want some uncertainty about the locations themselves? Anything (or everything!) that goes into the GP formula can be PyMC random variables, and you can have as many observed variables as you wish scattered throughout the model. So if you have a collection of noisy location observations for a single location over time, it would be easy to use that as observed data and estimate the underlying true value.