Turbulence data sets produced from computational fluid dynamics (CFD), especially from fine-resolved direct numerical simulations (DNS) and large eddy simulations (LES) of turbulent flows, tend to be very large due to high resolutions adopted to accurately resolve the smallest scales. While the computational capacity of high-performance computing (HPC) platforms has kept increasing, storage capacity has lagged to the point that more data is being produced than what can be efficiently managed. Among the several methods emerged to deal with this problem, an efficient technique is data compression. In this study, we present a proof of concept of a novel data compression approach that relies on Gaussian process regression (GPR) within a Bayesian framework to handle data sets in such a way that initially discarded information can be recovered a posteriori. The approach can be used to supplement existing compression algorithms with measures of uncertainty and we show that it can be applied to compress not only the 3D spatial fields of turbulence but also the discrete sets of time series data. The compression algorithm has been designed for data from spectral element method (SEM) simulations but can be extended to spatiotemporal fields obtained from other methods arising in engineering and physics. Our investigation shows that it is possible to use Gaussian process regression for data compression, however also highlights several of its limitations, in particular, that efficient implementations of GPR are crucial for its adoption, and that, while it is unlikely that the method can compete in terms of throughput with state of the art methods, given the cost of GPR, there is potential in terms of compression performance, as long as efficient bit-plane coding is integrated.