Abstract: In many applications, including climate-model emulation and calibration, there is a need to learn the conditional distribution of a high-dimensional spatial field given a covariate vector, based on a small number of training samples. We propose a nonparametric Bayesian method that decomposes this challenging conditional density estimation task into a large series of univariate autoregressions that we model using heteroskedastic Gaussian processes with carefully chosen prior parameterizations. We describe scalable variational inference based on stochastic gradient descent. The resulting generative model can be used to sample from the learned distribution or transform existing fields as a function of the covariate vector. We provide numerical illustrations and comparisons on simulated data and climate-model output.