  An important application of next-generation wide-field radio interferometers
is making high dynamic range maps of radio emission. Traditional deconvolution
methods like CLEAN can give poor recovery of diffuse structure, prompting the
development of wide-field alternatives like Direct Optimal Mapping and $m$-mode
analysis. In this paper, we propose an alternative Bayesian method to infer the
coefficients of a full-sky spherical harmonic basis for a drift-scan telescope
with potentially thousands of baselines. The can precisely encode the
uncertainties and correlations between the parameters used to build the
recovered image. We use Gaussian Constrained Realisations (GCR) to efficiently
draw samples of the spherical harmonic coefficients, despite the very large
parameter space and extensive sky-regions of missing data. Each GCR solution
provides a complete, statistically-consistent gap-free realisation of a
full-sky map conditioned on the available data, even when the interferometer's
field of view is small. Many realisations can be generated and used for further
analysis and robust propagation of statistical uncertainties. In this paper, we
present the mathematical formalism of the spherical harmonic GCR-method for
radio interferometers. We focus on the recovery of diffuse emission as a use
case, along with validation of the method against simulations with a known
diffuse emission component.
