Possible ways to speed up spatio-temporal GP modeling

The Kronecker stuff sounds like an option then!

The problem is that all observation locations of the satellite do not have data at a particular measurement time. In other words, depending on the measurement time, some locations have missing data.

How many are missing? If you use LatentKron (I think this will work for MarginalKron too but I’d have to try it) you can still pass a full grid of X_spatial and X_time. Then index the GP where you actually have observed data and feed that into the likelihood and then do NUTS sampling. You’ll get the result of the GP interpolating the missing spatial/time points “for free”.